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Chemical Phase Equilibrium and Physical Structure Within 
the Earth’s Mantle. 


By 


Yasuo SHIMAZU 


Institute of Earth Sciences, Nagoya University, Nagoya. 


Abstract 


Phase equilibrium of Fe-Si-O system within the earth is discussed. Due to the action 
of gravitational differentiation, there is chemical squeezing of the lighter material out of 


the earth. Then the silicate or oxide phase is decomposed below certain depths even if the 


temperature is low. It is concluded that FesSiO, layer is decomposed at 1,000 km from 
the earth’s surface. This result suggests the origin and the physical properties of Jeffrey’s 


20° discontinuity. 


$1. Introduction 

In my previous paper (Y. Sammazu [1955]), we 
studied the chemical equilibrium condition 
within the earth. Our discussions were based 
upon an assumption that a mixed crystal can 
be chemically stable with any mixing ratio of its 
components. For instance, FeO and SiO, were 
assumed to be able to combine with any mol 


in the combination, so that the mol ratio cannot 
diverge from some definite ratios. FeO and SiO, 
can exist either in the form of 2FeO-Si0O, (faya- 
lite)or in the form of FeO-SiO, (ferrosilite). Ina 
system of components which can combine into 
different compounds, gravity will obviously 
tend to develop a mechanically stable arrange- 
ment of different layers of different chemical 


ratio. There is, however, a very narrow area phases and of different densities. If the 


vapor 


water 


Fig. 1. Phase equilibrium of H,O-system in the gravitational field. 


gravitational differentiation proceeds far 
enough, as was discussed in my previous 
paper (Yasuo Sutmazu [1955]), then the system 
becomes saturated at the lower end with one 
compound phase, and further transport of 
material will form another compound. We 
then have two phases of different compound 
concentration ratios. There will be a concent- 
ration discontinuity at the point where the 
chemical potentials of the two phases are 


equal. 

It is interesting to note that ordinary Grsps 
phase rule cannot be applied to the phase 
equilibrium in the gravitational field as 
W. SwirtosiAwsky [1947] pointed out. He dis- 
cussed one-component system (H,O) under the 
action of gravitational force. According to 
the ordinary Gress phase rule, the maximum 
number of co-existing phases for H,O-system 
is three. In fact, at the so-called triple point 
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in the phase diagram we have solid ice, liquid 
water, and water vapor. 

Let us assume that large amounts of water, 
ice, and water vapor are put in a closed 
container. In Figure 1, the initial state (a) 
and final state (b) are schematically shown. 
First, the block of ice is partly in the vapor 
and partly immersed into the liquid. The 
immersed part is subjected to the action of 
hydrostatic pressure, which changes from the 
top vapor pressure f)(=0.006atm. in this 
system) to py+4p and thus the equilibrium 
temperature between the solid and liquid phases 
changes from T° (triple point=0.0078°C) to 
T°—AT as seen in the phase diagram. In 
Figure l(a), there is only one line along 
which T° can actually be found. As is easily 
understood, the ice below the surface of the 
water undergoes melting. On the other hand, 
the ice surrounded by vaporous space is 
covered by ice formed. Thus the movement 
of ice toward the surface takes place. The 
balance of heat requires that the heat used 
in melting the ice must be compensated by 
that in the formation of the corresponding 
amount of ice. As we see in Figure 1 (b), 
three phases have no mutual contact with 
each other in the final state. We may also 
note that the order of layering of the three 
phases is governed by their density differences. 


§2. Calculation of free energy of oxidation 


Along the general line in section 1, we shall 
discuss the chemical equilibrium of three- 
component system (Fe-Si-O) within the gravi- 
tational field of the earth. For this purpose 
we must get as a function of temperature the 
heat of formation 4H and Gress free energy 
change 4G during the chemical reaction. We 
can obtain 4H from the change of specific 
heat 4C, 


CT) 


Since the specific heat C, of each component 
is known experimentally as a function of 
temperature, 4C, can be obtained from the 
difference between the C, of each side of 
reaction equation. The procedure of calculation 


- > a reat 


is due to W. Erre. [1951]. Experimental 
values are obtained from the works of 
F. Bircr et. al. [1942]. We calculated 4G° (=4G 
at ordinary pressure) from the integration of 


sale (2) 
Ol No: T? 
Results of calculations are as follows 
Fe+10.2FeO | 
4H=—65.6 x 10°+0.631 T+0.550 x 10-* T? 
(cal/mol). 
4G® = —65.6 x10°+23.2 T—0.631 TinT 
—().550 x10 7? (cal/mol). (32K) 
Si+0.2Si0, 


4H=—221.9 x 10°+8.79 T—0.26 x 10-3 T?, 
AG® = —221.9 x 10?+116-5 T—8.79 T inT 
+0.26 x 10-3 T?. 
2FeO+Si0.Fe.SiO, 
4H=—343.0 x108—2.45 T+1.76 x 10-3 T?, 
4G = —343.0 x 10?+64.9 T+2.45 TinT 


(4) 


—1.76 x10-8T?. (5) 
Table I. Heat of formation 4H at several 
temperatures (kcal/mol) 
Reaction Reaction Reaction 
TC (3) (4) (5) 
1,000 —64.42 ~213.4 ~ 343.7 
2,000 — 62.14 — 205.4 — 340.9 
3,000 —58.76 —197.9 — 334.5 
4,000 — 54.28 — 190.9 — 324.6 
5,000 —48.70 —184.5 —311,.3 aye 
Table I]. Free energy change 4G® at 
several temperatures (kcal/mol) 
Reaction Reaction Reaction 
TCR) (3) (4) (5) 
1,000 —47.3 —165.9 262.9. 
2,000 —31.0 —121.5 — 183.0 
3,000 — 16.2 — 82.6 —105.8 
4,000 — 2.53 — 43.4 — 30.3 
5,000 + 9.78 =e AT + 41.8 


Since no experimental 
another reaction 


FeO0+Si0,2F eSiO; 
is not discussed here. 


value is available, 


(6) 
Tables 1 and 2 give 
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4H and 4G® at several temperatures. The 
quantity 4G° is nothing but the chemical 
potential at the standard pressure. We have 
also the following general relation between 
chemical potential 4G and the vapor pressure 

4G=4G°+RT InP. oe 

The equilibrium condition is given by 4G=0. 
We get the temperatures T°=4,200°K, 5,200°K 
and 4,400°K at which the equilibrium condi- 
tion are realized at room for the reactions (3) 
(4) and (5), respectively. At these tempera- 
tures the corresponding oxides or silicates are 
., decomposed. At higher temperatures, oxides 
or silicates are unstable. 

Pressure dependence of the decomposing 
temperature J° is given by the CLapEyRON 
equation 

dT’ T,4V 
ae AR 
where 4V is the volume change between the 
left-hand and right-hand sides of the reaction 
equation. In (3) and (4), 4V are equal to the 
volume changes in gaseous oxygen phases, or 
—3X22.4 x 10’cc/mol and —22.4x10%cc/mol, 
respectively. In (5), 4V is the difference in 
molar volume between both sides of the 
equation, viz. 
2FeO+Si0, @ Fe.SiO,, 


(8) 


4AV=—3.62 cc/mol. 


23.96 27.27 47.61 
In conclusion, we get 
at? 
{OG dP =1.76 deg/atm. (9) 
for (4) =1,36 (10) 
ToL Ca) =1 19x 10-°. (11) 


These equations show that an exceedingly 
high temperature is necessary to decompose 
oxides within the earth’s interior. Then we 
may infer that oxides and silicates are stable 
within the earth’s mantle if the abundance of 
elements concerned is appropriate. 

§3. Phase equilibrium within the gravita- 
tional field. 

As we see in §1, the phase equilibrium 
within the gravitational field has a particular 
character, that is, the ordinary Grsss phase 
rule cannot be applied. In one-component 
system, as stated in §1, the movement or 


migration of solid ice occurs. The melting of 
ice at the bottom has a compensated formation 
at the top. In a multi-component system, on 
the other hand, there must be the migration 
of any one phase through others. This is 
nothing but the diffusion phenomena. Even 
in crystalline solid system, there will always 
be some diffusion effect. Thus, the chemical 
potential gradient will always cause a diffusion 
through the system if sufficient time is allowed. 
We do not concern ourselves with the mecha- 
nism of diffusion. The final equilibrium state 
is assumed apviorz. 

Migration is controlled by gravitational 
differentiation in which the lighter material 
moves upward and the relatively heavier 
material moves downward. If the lighter 
material is in the fluid form, diffusion through 
the solid phase will effectively realized. Equi- 
librium of hydrous minerals, such as mica, 
has already been discussed by several authors 
based on this concept of migration (cf. J. B. 
Tuompson [1955]). Of these studies the most 
famous is that on calcite (CaCO,)—wollastonite 
(CaSiO3) reaction. (cf. ‘‘Theoretical Petrology’”’ 
by T.F.W. Barra [1948] (p. 286)). As this 
reaction proceeds to the right, gas phase CO, 
is produced and the equilibrium -temperature 
will therefore increase with increasing pressure. 
This is given by equation (8). At one atmo- 
spheric pressure, the formation of wollastonite 
requires a temperature of about 450°C. 
However, the rock is usually semipermeable 
to CO, gas which thus escaping from the 
locale of the reaction will percolate to the 
surface. If the gas escapes as soon as it is 
produced, the molecular volume of the gas can 
be neglected and the change in the temperature 
of the reaction with pressure can be calculated 
from (8), taking into consideration only the 


molar volumes of the solid phases. The 
molar volumes in c.c. are as follows: 
CaCO;+Si0@CaSiO;+CO: , 

AV=—18c.c./mol , (12) 

di” = mh /208 —(0.013 deg/atm, (13) 


dP 106 


the latter showing that the differential pressure 
squeezes CO, gas out of calcite at a lower 


~ se 
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temperature than is otherwise possible. The 
formation of wollastonite is facilitated with 
increasing pressure. In this case gas escapes 
through the fissure rather than by diffusion. 
As is shown later, however, the above discus- 
sion is not rigorous since the chemical poten- 
tial gradient of CO, is neglected. 

The idea ‘of chemical squeezing of some 
light materials out of deeper parts of the 
earth has also been presented by several 
authors. Decrease in the degree of oxidation 
or decrease in the ratio of oxygen to metallic 
ions at different levels in the earth’s crust 
shows the squeezing of oxygen. H. Brown 
and C. Parrerson [1948] conclude that the 
main variable elements in meteorites are iron 
and oxygen, with the decrease in oxygen al- 
most exactly matching the increase in iron. 
It proves that oxygen is squeezed out of the 
deeper parts of the protoplanet of meteorites. 
T.F.W. Barta [1948] calculates the phase 
equilibrium of Fe-O system in the crust and 
concludes that magnetite (Fe;0;) is unstable 
in the deeper parts and underlain by a layer 
of wiistite (FeO). At still greater depths, Fe 
becomes stable. 

As our first step, we shall study the general 
theory of equilibrium of reaction 


C 


B+G 


Fig. 2. Phase equilibrium of the reac- 
tion A@B+C in the gravitational 
field (08 >04>00). 


A@B+CT, (14) 


where A,B,C, for example, correspond to 
FeO, Fe, O, in (3). We assume that C can 
migrate through the system. At low pressure, 


the reaction (14) proceeds to the right if the temp- 


erature is higher than T°. If os>o4> 0, the 


thick layer of A is unstable in the gravitational 
field even if the temperature is below T° 
in (8). A more stable state exists if the 
reaction proceeds to the right below a certain 
depth h, so that the dense phase B is formed, | 
while the liberated phase C is displaced | 
to the top of layer A. The top of layer A 
is the zero level surface. This state is shown 
schematically in Figure 2. The migration of 

C phase is shown by a pipe line (fissure) 
which passes through A layer. 

Gisps free energy G per mol is a function 
of pressure P, temperature T, and the depth 
from the surface h. Thermodynamical iden-” 
tities give 


OGi er 
(ar aoe oe Se 
OG: a 
Gaye coe! 
OG: 
—— =—Mg, 
(age) g ee 
where Si, Vi,M:i are the entropy, volume, 


molecular weight per mol, respectively, and 
suffix z denotes A or B or C. The assump- 
tion of hydrostatic equilibrium gives 


ibe mel 
ays =p4g— Vi g 


(18) 


in A layer. Incompressibility is assumed which 
is allowed in our chemical problem. As we 
see in the theory of chemical equilibrium, A 
phase, for instance, does not mean that there 
is no B or C phase at all. It means only ™ 
that the partial pressure of A is predominant. 
Though B phase at the depth hk is under the 
pressures P=p.gh, this is not so in C phase at 


depth h. The pressure P’ of C is given by 
Che M, 
an °9= 7 9 (19) 


Schematically P’ is the pressure in the fissure 
in Figure 2. Thus the phase equilibrium at 
depth hf is attained even if PP’. This 
circumstance is similar to osmotic equilibrium. 
The wall of the fissure corresponds to a 
semipermeable membrane. Then from (15)-(17) 
we get 


dG 4,rn= —SidT+ VadP= —S,dT+ Vuoagdh 5 
dGa,n= —SszdT+ VedP= —SzdT+ Vzosgdh ; 
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Ge n= —S,dT+ V.dP' = —S,dT+Megdh . (20) in which Ss+Sc—S,4 equals (4H/T) by defini- 


Chemical equilibrium gives ~ 


AG 4,n=dGxn+dGe,n . CAL) 
Substituting (21) into (20), we have 
aT 9 —Va0s+Vu04+M.) (22) 


dh Sz+S.—Sa ; 


adiabatic! 
1} 


\ 
3 
Fig. 3. 
(a) 


FeO>Fe+30, 

(b) SiO,>Si+Os, 

(c) Fe:Si04>2FeO+SiO. 
adiabatic: 
melting: 


is obviously zero. Even if the temperature is 
~ below JT, for the whole system, A phase is 
decomposed into B and Cat depth h. The lower 
the temperature, the greater the thickness h. 
We shall now make the calculation for 
Fe-O system in (3). As o4=6.00, 02 =7.87, 
Mzg=55.8 and g=10? in the mantle of the 
earth, we get Curve (a) in Figure 3. Similar 
calculations for Si-O, and FeO-SiO, systems 
are also given in Curves (b) and (c) in Figure 
3, respectively. We see that SiO, is very 
stable and will not be decomposed within the 


earth. 


tion. Since Mi=Mz+M, we get finally 


ap TMug(*—1) 

dh 4H : 
We see that the thickness of A layer is 
determined by (o4/0z). At temperature T°, h 


(23) 


sol 


XHO™ kK 


melting 


\ 
\ 


Temperature-depth curves of chemical stability. 


temperature derived from adiabatic gradient 
temperature derived from melting point gradient 


If the condition px >e4>pc is not fulfilled, 
the state given in this section does not occur. 
The equation (23) can easily be modified 
T(4V.+ £4Ve) 

OA 
4H 


Lie 


aie 4 


where JV; is the volume change of the solid 
phase (=Vz—Va) and 4Ve¢ is of the gaseous 
phase (V,). Barta’s discussion in (13) cor- 
responds to the approximation in which the 
effective volume change in the movable phase 
(oc/pa)dVe is neglected. The equation (24) is 
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a generalized CLAPEYRON’s equation corres- 
ponding to (8). 

Though the discussion in this paper does 
not give the sufficient condition for the equili- 
brium, it does give necessary condition. If 
the abundance ratio of each component 
is not suitable, the thick layer A cannot 
be formed and thus our conclusion fails. 
From cosmic, meteoritic or petrological 
abundances given by several authors, we may 
safely assume O: Fe: Si=49: 34: 17 for 
the simplified model of the earth. Since 
SiO. is much stable than FeO in the earth 
SiO, is first formed and the remaining oxygen 
combines with Fe to form FeO. The above 
abundance ratio is sufficient to form a thick 
SiO,-FeO layer. The excess abundance of Fe 
will make the metal core. Ata certain depth 
Si and O, will be squeezed out of FeSiOs 
lattice and Fe.SiO,(fayalite) lattice will be form- 
ed. The critical depth for FeO-stablity (Curve 
(a) in Figure 3) is much shallower than that 
for Fe,SiO,-stability (Curve (c)). Then, silicate 
mantle (FesSiO,) is split into Fe region at a 
depth greater than Curve (c). In Figure 3, 
the temperature-depth curve for the mantle is 
also given. ‘‘ Adiabatic’’ is derived from the 
adiabatic temperature gradient and “‘ melting ”’ 
from the melting point gradient. The former 
gives the minimum plausible temperature and 
the latter the maximum within the earth. (cf. 
Y. Suimazu [1954]). FeO,SiO, and Fe,SiO, 
are, thus, stable at the region shallower than 
the points at which the temperature-depth 
curves intercept curves (a), (b), (c) respectively. 
The instability of FeO layer at some depth 
has already been shown in my previous paper 
(Y. Samazu [1955]) from a different point of 
view. 

If the earth is formed by accretion process, 
the oxidized or silicate compound will be 
stable in the earth’s initial stage. As the 
accretion process proceeds, the silicate phase 
sphere grows and the temperature is raised at 
the same time. Then the critical depth as 


shown in Figure 3 is reached at some stage 
and a discontinuity in composition will be 
formed. Possible existence of Fe region due 
to chemical squeezing at depths greater than 
1,000 km suggests the jump in electrical con- 
ductivity there and the existence of JEFFREY’S 
20° discontinuity. 

This work was carried out by the writer 
under receipt of Grant in Aid for Scientific 
Research of the Ministry of Education. 
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On the Variation in Bulk Modulus/Density in the Mantle. 


By 


Michiyasu Saima 


Geophysical Institute, Kyoto University. 


Abstract 


Qn the assumption that the mantle is composed of the ionic crystals, the distribution 
of the ratio of bulk modulus F to density o and of pressure P is investigated by means of 
atomic theory to reveal the character of change occuring between the layers B (33-413 km) 


and D (1000-2898 km). 


It is found that, if this change were the polymorphic transition 


from the low pressure phase to the high pressure phase, the decrease of the gradient of 


change may be of chemical composition. 


k/o would occur between B and D. This is against the observed result. 
If so, the dissociation energy MH and the reduced 


Then, the 


density 9 increase and the inverse power 2 of potential between the atoms hardly changes 


there. 


§1. k;/o derived from the seismic P and S 
wave velocities within the earth increases 
steadily in the so-called B- and D-layer. Then, 
in studying the density distribution BuLLEN 
(1947) assumed that the materials are of 
homogeneous constitutions in B and D. In 
this paper, this assumption will be adopted. 
But s/o increases abnormally in the inter- 
mediate layer (C-layer) between B and D, 
and such an abnormality has been supposed to be 
due to changes of constitution of the materials 
in this layer. For example, 1) Bernat (1936) 
suggested that such changes could be attri- 
buted to a rearrangement in the crystal 
state of olivine to the high pressure phase 

and 2) Birca (1952) proposed that it could 
be attibuted to changes of chemical compo- 
sitions. Here by means of the theory of matter, 
we shall firstly study the possibility of the 
polymorphic transition. 

§2. Since the materials which constitute 
the earth interior are composed of many 
kinds of minerals of complex crystal structures, 
it is difficult to treat the materials as such. 
Then we are obliged to reduce them to the 
simple ionic crystals representing their 
characteristic properties. In the following 
study, NaCl and CsCl types of crystals will 
be taken as example. NaCl type crystal changes 
polymorphically to CsCl type at more than 
several thousands atm.. For example, the 


polymorphic transitions occur in AgI at 4500 
atm. and in Nal at 50000 atm. (Jacogs, 1938). 


O a 

(009) (200) 
CsCl type 

Fig. 1. Crystal lattice 


As the modulus of elasticity obtained theo- 
retically as a function of pressure and tem- 
isothermal, the adiabatic ks/o 

the seismic waves velocities 


perature are 
derived from 
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must be converted to compare the calculated force and the repulsive force due to the | 


elastic modulus with the observed one. The overlapping of the closed shells are assumed be- 

tween the different kinds of ions, and the 
ks/kr=14+ Tare , (1) repulsive electrostatic force of the above kind 
between the same kinds of ions. Then we 


ratio is 


where ks is the adiabatic bulk modulus, kr the 
isothermal one, I the temperature, a@ the heve 


coefficient of thermal expansion, y¢ GRUNEISEN he Os by pena t Ra) 
constant. With depth h, T increases and @ YAB UL weh™ 

and y¢ decreases slowly, whereas Tay¢ seems WARE 

to be constant. Assuming T=1000°K, a=2.0- bav=( Za 402 +1) (4) 


10-®, ye=1.5, we have kr which is smaller 
about 3 per cent than Rs. 
kr/p and p are expressed as follows (Suma, 


where Nu, Nz the numbers of outer shell 
electrons, Ru, Rez their ionic radii, Z4, Zz the 


1954); the second repulsive term derived by PAuLING 
kr/o=len+2Ci2+P)\/90 ; ey (1940), he reduced the ion to the hydrogenlike 

Pp ge Dea)e=o ; one and simplified the expression of potential 

vi ae derived by the MHerrter-Lonpon method, 

F=0+- N > loe( 47") i (3) assuming that the ions have the certain 

tin Aaa radii Ry and Rz, and the repulsive potential 


where JN is the numbers of atoms per mole, . : : 
increases, when the distance between the ions 
k BoutTzMAnNn constant, v the volume per atom, = 
; -? becomes less than the sum of these radii. 
wo, frequency. The attractive electrostatic Thea i 


(aes ar ote | Fess B Ty), 
2uN 7 Te i fie ZINE (oh 


1 = AL? + 1? + 13”) + @{2(132e1 + 1972+ 132€3 + Lol sy + Iglies + Llees) + 12(€1?+ €5?/4 + 62/4) 
t+ 127(€2? + €g7/4 + €4?/4) + 132(€3? + €52/4 + €42/4) + (Jol3/4)(€s€5 + 2€2€1+ 2€,€3) 
+ (131,/4)(€s€5 + 2€1€5 + 2€3€5) + (L,12/4)(€1€s-+ 2€1€5 + 2€1€,)} 
++ +- ++ +- 
BT»° = (ee + +B_-)S,+28, eS == B(Sp+ OSn) = (ge ; ( 5 ) 


where /,,/,, 2; are the coordinates of the lattice points taking a in Fig. 1 as unit. Using the 
potential, we have the vibrational energy terms as follows (BRADBURN, 1943, Gow, 1944); 


kT 3x hos WORT 1 oe h 
S16 Cer) = (Hel dats] 2 
g og det [27] 9 og u+log rT)? 


vN p=1 kT 6 
| dir! pr" $3" 
det [¢j7]=\ 21" doz’ do" 
/Par' baa’ 3s" |, 
1 ag? 10/0 49 | 
ji O47 i ae 
: : NT On ade Ti alae a) (oe 
where / is the mass of an atom. Then from (3), (4), (5), (6), p and ci; are 
RG Cc p \&%-1/8 (R3?+2R3'+5R,!)T 
p=— of (S>—r( ) a) ) Eo SURES E 
A DCS 00 nea) 2(Ri°+ R38 : ay 
RO = (n-1)/3 
(Cim= A als {S.—S, +7 mh (Toei Ieee) | 
{(Ri + Ry')(Ry?+2R3"+5R,')/2+(R!+R)(RE+2R,+3R32)/8 
+ ss. PRN RY+2R5")-+(7(Ra' P+ 2R2Rs4 + (Rs)?+- 2Ro!R32)/4}T |. (8) 


2(R+R,)20 


ionic charges, B)a constant. With regard ta, 
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RO n= 
C2 ries (s. asl oa) a 2 | 


fmere, NRE RIM SRI 2+ (RL IR NY +2RPR3!/4 
- «abel ose leniniiag _+2R3'(Rs?-+2Rs") +R YRS 4 9(Ry! 1, 
(Ry +R )0 ; stay a 
= kb 0 \rvenn pm hee 
n—1\ 0 ia a ke 
1,2? 2q]..27r ©, a Siacd 
s,nunas  hirhtls (1,0, 0) (0, 1, 0) (0,0, 1) > (1) 
(Ly? + Da? + 1g2)"/2? (2,0, 0) (0, 2,0) (0,0, 2) > (2) (10) 
CL, 7; OF (5-0) 1-0, 15 1) 2X, 1, 
Rerr=[(—18-1-1+2)---.426—p) 9p, +7( 2)" 
Po 
(—It-1-(42)-- fn. 2(s—)} Typ, | BE (2 aces 
a Po 
Sn=3S),, 352,+6S92,=S,,  S$@,429¢2 =9@, 25054 sam— sa» (11) 
N+ n+ N+45 


where , is o at absolute zero temperature and zero atm., RO=E the dissociation energy, R 
the gas constant, A the mean atomic weight. S,@” is represented by integral as follows 


. 
S,0M = ~ ur/2-1 1,2],24 2rp-l2u = 1 
Pinay RPMI dam 
The fundamental functions o and o®@ are represented by @ functions as follows 
OQo1 4  (4u) = 833(0, e~™) + 363(0, e-*).2(0, e-) —1 
Tn 4 (44) = 823(0, e-) +3037(0, e-™)02(0, e-*) , 

0705(0, e-4#) 
Ou 2 
070,(0, e~* 
Ou 
oP. = 9590, 2-8) —1, o{22,,, (dt) = 0.9(0, €-*) , 


—4U 
= 95%(0, ey). | ae, =82%(0, 68") 


[i yr To@")(n)du . (12) 


0 aoe Be 00. Se es 


a2, , (41) =05?(0, e-*) +6. (0, e-*) ") 4.20,(0, €-*)85 (0, e-4%)° 
0? oth ent 


Ou? 


nis ee) ; 


)420.(0, e-*), (0, e-!*) . 
Ou? 


a Poi» (400) = 0,0, e-*) 4.4.20, e-*) 


pO a + e % « 7 =) 
; Ou? 


(13) 


where we use the Born’s method (1940) in the calculation of S,“’"”. Since, firstly, the pos- 
sibility of the polymorphic transition from the lower pressure phase to the higher one is in- 
vestigated, NaCl type crystals of the lower pressure phase will be taken as the materials in 
the B-layer. Then, adopting the BuLien’s density distribution in this layer, we can determine 
the constants in the expression of kr/p and p corresponding to these materials. The adiabatic 
temperature distribution derived by VERHOOGEN (1952) and the mean atomic weight of 
magnesium olivine (A=21) is assumed in this layer. Then 
o=1, 20, 8= 225000, po=a-45. 

In this transition, the atomic configuration changes, but the constants 7, 6, 6, 09 do not. The 

density incresaes from p to 0’ by about 4 per cent at the polymorphic transition, and we have 


1/ 
Pr_-2.969-100| {—1.0213( ° y "+2.298(° y+5.333-104 71 h 
Po Po J 387 


0 
/\1/3 / pao 
+{—o.9440(-" ) < 1,989( p ie 5.333-10-5 TL” a . (14) 
{ ry) 09 ) 387 
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If the transition occurs at the same time 
between 413 and 800km, k,/p varies as shown 
by the dotted line in Fig. 2. We see in Fig. 2. 


—— Observed 
teccee: Calculated 


1000 


k7/o in the mantle 


1500 = (km) 
Fig. 2. 


that the calculated and observed distribution 
kr/o have an opposite tendency, and thus k7/o 
derived from the observed value denies the 
possibility of the polymorphic transition which 
was suggested by Bernat (1936). Next, the 
possibility of chemical changes will be studied. 
Only the changes in the constitutive atoms 
of NaCl type ionic crystals, for instance 
CaO — MgO, will be investigated here. Since 
kr/p is smooth in the boundary between C- 
and D-layer, the materials in the lower part 
of C-layer seem to be the same as in D-layer, 
and the value of the constants may be con- 
tinuous. Considering these conditions we 
carry out the above process again and get 


d=1, 0@=225000, po=3.45, m=6, A=21. 
6=1, 60=430000, po=4.20, m=5, A=21. 


Namely, the dissociation energy & and the 
reduced density go, remarkably increase be- 
tween 413 and 1000km, whereas » hardly 
changes there. The density distribution is 
determined so as to agree with the observed 
total mass and the total moment of inertia of 
the earth and the error of the estimated 
density at 1000km is less than 0.1. With 
regard to the effect of the density distribution 
between 413 and 1000 km on the pressure at 
1000 km, its deviation from BULLENN’s pressure_ 
is less than 2 per cent, assuming a linear 
density distribution instead of BuLLEN’s. Then 
his density may be adopted to determine the 
crystal constants. If the temperature gradient 
is steeper than the adiabatic one, the density 
gradient is smaller than that of BuLLEN. But 
as the term dependent on the temperature 
cancels the term independent of it, this 
influence is negligible. 
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The Dynamo Theory of the Earth’s Main Magnetic Field 


By 
Hitoshi Takeucut 
Geophysical Institute, Faculty of Science, Tokyo University, Tokyo. 


Abstract 


The possibility of the dynamo action in the earth core is considered in the case when 


the toroidal 


fluid velocity is larger than the poloidal velocity. 


Infinite numbers of 


stationary dynamos, each corresponding to a different eigen-value y=4rneaV»..,, are shown 


to be possible there. 


The dynamo including the dipole field is that for the smallest y. 
The geophysical significance of this result is considered in section 4. 


A mechanism of rever- 


sion of the earth’s main magnetic field is proposed in section 4. 


§1. It is now suspected that the earth’s main 
magnetic field is maintained by the dynamo 
action in the earth core. The essence of this 
dynamo theory of the earth’s main magnetic 
field is as follows. There may exist in the 
earth’s core chains of relations in which a 
field H, interacts with a motion V; to produce 
a field Az, which, in turn, interacts with 
another motion V, to reproduce Ay. If it is 
possible to show the existence of self-sustain- 
ing dynamo in the earth’s core, the main 
field may be explained. All that is necessary 
is a weak initial field to start, and the field 
will then grow automatically until the_energy 
dissipated becomes equal to that supplied. 
Thus, the central problem here is to examine 
whether Maxwe.u’s equations of  electro- 
“| dynamics possess solutions representing such 
a self-sustaining dynamo. 


A method has recently been developed 
(Exsasser 1946, 1947; Taxrucut and SHiImazu 
1952, 1953; Bunuarp and Gentuman 1954) for 
deciding whether any specified motion in a 
fluid sphere will act as a dynamo. In this 
method, solutions are obtained by expanding 
the velocity and the fields in spherical har- 
monics to give a set of simultaneous linear dif- 
ferential equations among the radial functions 
corresponding to the respective fields. As the 
plausible fluid motions in the core, two types 
of motions are selected. One is of rotational 
type and the other is of convective type. 
They are usually called the 7,° and S,” 
motions, respectively, and are schematically 
shown in Fig. 1. Mathematical expressions for 
them will be obtained by putting »=1, m=0, 
T(n)=V%(r) and n=2, m=2, S.2(r)=V2?(7) 
respectively in the following expressions for 


S,o T.0 Tze Vo V2 
N N N N $= 
Ro a 
‘ae . pe IG 
s s S 5 ¢=0 


Figure 1. 


the (7, 8, ¢) components of poloidal and toroid- 
al type velocity or fields. 
S type (poloidal type) motion or field 

y component 


= —n(n-+ Sone Yn™ , 
6 component 


Se m |an—l OYn™ 
=-| op EGAN } AYA ; 
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o@ component 
0 Yn 


aSnt ( m pa ee 
=-|7 dr +(n+1)Sn |r sin 0 0¢ 


T type (toroidal type) motion or field 
y component=0 , 


, (1.1) 


ONE 
6 component=—Tn™(7)7” sin 008” 


Mantle 
Figure 2. 


OYn™ (1.2) 


@. component= Tn™(”)r”— a0.” 
where 
Yn™=Pn(cos O)gin M0 - (1.3) 


The upper suffix 2c in S,.?* means cos 2¢ in 
Ye 
In order to study the maintenance of magne- 


tic fields by electro-magnetic induction, we 
must solve the following equations 
divH=0, (1.4) 


Figure 3. 


OH 


4rkK By a eae COV 8 B= 0, (1.5) 


0 
where H and V are magnetic field and fluid 
velocity. The system for which the solution 
is required is shown in Fig. 2. The coordinate 
system is fixed to the uniformly rotating 
mantle which is assumed to be an_ insulator. 
Thus V in (1.5) denotes the fluid velocity 
relative to the rotating mantle. By the 
assumption of ‘‘viscous’’ fluid, V;,V, and Vy 
must be zero at the core boundary 7=a. 

Taking (1.1) and (1.2) into consideration, 


we can satisfy these boundary conditions by 
putting 


dV..?(r) 
+ =() 1.6 
soa), ae) 
at y=a. In the numerical calculation to fol- 
low, we shall take rather tentatively 


Vir 1r) = V22"(7) — 


Vin) = ; VioViE), V22%(r)= <Va.nVi( 


Vi=1-€,. Vir Saal ee 4a), 


(1.7) 


which satisfy (1.6). In (1.7), Vi.9 and Vane 


_.- respectively. 
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are some undetermined constants. 

On assuming the existence of the above 
fluid motions in the earth core, we have there 
a group of magnetic fields which are mutually 
coupled by the inductive action of the fluid 
and among which the dipole type field S,° is 
included. The way of coupling of the fields 
with 2, m<4 is shown in Fig. 3. In this 
figure, the coupling between rows is by the 
S.°* motion, that along rows by the 7," motion. 

It is needless to say that the chain is 
continued infinitely with infinite numbers of 


magnetic fields. In Fig. 3, S,° and T,", for 
example, mean the magnetic fields by putting 
m=1, m=0 in (1.1) and m=2, m=0 in (1.2), 
respectively. The first four fields in the chain 
are shown schematically in Fig. 1. 

§ 2. It can be shown (Takeucni and Suimazu, 
and Buntiarp and Getuman, loc. cit.) that the 
whole problem is reduced to solve the simul- 
taneous differential equations of the following 
types among the radial functions correspondig 
to the respective fields 


d*5," doe = 9x 74 sy [28 
E dé? +4 dé yEViT, 

Gai dT,° Se 2 dV,° C 0 72 ¥ 9 2 9 iq . 
bee oar eet ae tS alee Vite Ty?) + gee VIELE | 

aT, dT 2 6 
fe ca 7 pg Vers ae ie edie |- bins 33 

@T25 —dT28 2 y d( dS," 

: 95-7. Te fe 3/7, 2¢ 1 317.20 

oo ae droke ora seek Vitae get tt af ey, 1s | 


under the boundary conditions 


T2°(€)= 1¢%])= f BE exh 0, 
0 gE ds,° 3 
S; @)+, dé =A (5.9) 


at f=10) in (2.1), SnXe), Ta), v2’ and y are 
Sr™E)=a"*Sn™(7), Tn E)=a"Tr™(r), 


z=4rKraV,.o, T= ATE aeons 


(2:3) 
In (2.1), only the magnetic fields 
with #<2 are taken into account. This seems 
to be a kind of theoretical cut off. Some 
reasons, however, will be given in section 4 
for that this is not a theoretical cut-off and 
by solving (2.1) we are studying the group of 
functions corresponding to the first eigen-value. 
From the above equations, we see at once 
that our problem is an eigen-value problem 
for « and y. The physical reason why we 
get the eigen-value problem here is as follows. 
From (1.5), we can easily get 


1 | rav=\v AHxXV)dU, (24) 
K 


where dU is an element of volume and _ the 
integral extends over the whole volume con- 
cerned. The left and right hand sides of 


(2.1) 


(2.4) denote the JouLE heat and work done on 
the field by the fluid, respectively. Thus, in 
order to sustain the fields against the JouLE 
dissipation, we must have a certain amount 
of work done by the fluid, and this certain 
amount is determined by solving the eigen- 
value problem for w and y. 

The eigen-value problem (2.1) with (1.7) and 
(2.2) has been solved numerically by TakEucal 
and SHimazu, and ButiARpD and GELLMAN (loc. 
cit.). The noted results obtained are as 
follows; 

(a) We get real eigen-value y in the whole 
real range of the ratio 


ey say. (2.5) 


Thus, the self-sustaining dynamo is possible, 
at any rate up to the point where spherical 
harmonics of the first and second degree are 
included in the argument. 

(b) Except for small a, the eigen-value ¥ 
and the corresponding eigen-functions vary 
very slowly with a. 


(c)’ When xz or a>, putting 
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14 
T?@) =a (o) (2.6) 
we have 
Sea VT, | 6 dT, 2 dVi'g, 
T; OS la rt dé? tie 4 dé ~ 3E dé re ’ 
dS," 4 dS _ 2724 pop, 3 “I y,2il? (2 dV.2 “yaver) Ts], 2.7) 
de ie ae Pas pe EV, dé == cae 2 


Putting (1.7) into (2.7), we get 


aT,” 12 dT. Ree aT,” 
Ge Teac wae Tg: dé 
4 6X ya pre| Fa é) a ie. 116)T |=0. 2.8) 
The corresponding boundary conditions are 
T dS? 3 go OTe 119 PTY 5 1g ATs" _ 9 2.9 
T.(E)=0 , dE Sa ae dé? + de? dé (2.9) 


at €=1. The eigen-value y in this case is 
about 42=y2, say. The eigen functions cor- 
responding to this eigen value are shown in 
Table 1. 


Table I. 

i é S19 T.9 T 28 
0 0.2019 — — 
On 0.2153 0.01524 1.8432 
0.2 0.3442 0.01254 1.1541 
O73 0.6316 0.028203 0.3900 
0.4 ().8477 0.021227 — 0.3454 
0.5 ().6290 — (0.027368 — 0.6733 
0.6 —0.03436 — 0.01427 — (0.3616 
O77 — 0.6540 — 0.01626 0.2767 
0.8 — (0.8300 —0,01301 0.6672 
0.9 — (0.6802 — 0.026757 0.5410 
1.0 —0.5 0.0 0.0 

It must be admitted here that the values 


of 6 &€7,?' (€) in Table V of our paper (ID 
were in error and must be divided by 6.99860 
to get correct values. 

§3. The results in the last section show 
that the self-exciting dynamo is possible, at 


any rate up to the point where spherical 
harmonics of the first and second degree are 
included in the argument. As was shown in 
Fig. 3, however, there are infinite numbers of 
magnetic fields with », m>2 belonging to 
the group we are considering and it may occur 
that the self-sustaining dynamo is impossible 
when the fields of relatively higher harmonics 
are taken into account. In order to make 
this point clear, we shall study whether the 
self-sustaining dynamo is possible or not when 
spherical harmonics up to the degree four are 
included in the discussion. In this case, we 
have twelve equations of the types as in (2.1) 
among twelve radial functions S,°(&), T.°(&), 
Ls?°(E), TPE), SaXE), Sst(E), Ss bya” TE), 
TPE), TAE), TYE) and T(E). When x or a 
tend to infinity, putting 


T(6)=2T YO), TP(E)=cT LE), 
we can transform these equations into 


Oe See Thee T,*¢ and T,#~(1/x)=0, 


(3.1) 
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8 /des® 4 dS,° 96 x9 we 
A ae te de )= 3 ITE 
Qe 

st 96 x 12 7 ie 


7 dé 


a ( @T! 6 AT,’), 48x6 
5 Ee Ras Fe a) 35 


288 x6 _ dV.2¢ 
T,? 2 720 
35 | ‘ (¢ a ) ae g 


192 x6 Cee 32 50 = dV.2¢ 
G20) Sat bey ee 2 0 UVa c \Ee 
fame Oe eek ta Ry Wf 7 (¢ aes eee Nee 


dV, 
S3° y => 
“ dé é 0 ’ 


° 


<b 
dé 


288 x 6 
35 7 
+ asi(E ave 


dé 
7x 6 yT2° Vee = 


eT: 7 ode 


aes ve) | es Bee ues V.2¢E? =(0 § 


de* 


x 12 ¢ 0 yo— 0x8. yT V2°E2=0 , 


’ 


16 x 20 dV,° 1. 80/d?T,°. 10 dT.2 
a = Nee AS) = — Ma bi pe ea 
eee Tae pA ee >) 
960, j= dV22¢ 4 dia 
cathe pga lig pa a0 7 2C OV7, 2064) _ 
: | “(E ie t3Y )43 9 ees 


480 x6 7 dV,?¢ : Close 1 
ee aA 10 | 0 2 V,2¢ 2 0 V.2¢ 6 


3 aS 


‘i 


160 x 207° V°E7=0 - 


Using the third and fifth differential equations in (3.2), we can express 


oy oye on yee 
ge TAS Vers, Cece +3V. ) Le 


ST vive 4 i a aT yee 


d (TY V22¢E8) 


pees ss 


15 


0, 


sh. 


Yd gt igre =a EEL a 8 ass? \+ 288 x2 y |-3 ( ist dS; At. is \V 


(3.2) 


Iho evaal See! shat 


terms of 7,°(&) and 7T,(&). Inserting the results thus obtained into the first differential 


equation, we get 
aS YE) , 4 dSi° 


eee TEU ek 


S,°(&), which satisfies (3.3) and which is finite at £=0, must be constant. 


condition for S,°(&) is by (2.2) 
GS ose 
fe AN) 
eee 
at E=1. From this, we have 
De) = const =0 © 


From (3.2), we get also 
d°S;? 8 dS; _28,,Vs" 
dau Fe devi Vio 
Gy MeO bee 1a led VY" 


= Ge ()) 
(fe fede 7 Eade 


(3.3) 


The boundary 


(3.4) 


G3) 


9 dT, ¢ on d V.7° Tr 
2 y2e V2" Tee =O 
r| SEV. dé +(33 ieee ae ) wy 


(3.6) 


When V,°=1—€ and V.?°=(1—€)? as in (1.7), the above equations can be reduced to 
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aT 2007 88 d*7,° 26/1 _£)2 [3 i dTe x ] 
: +(33—49&) To |=0 . oo, 
ae te ae te ae += “7 =F) Sel=6) dé ( ayy: (3.7) 
The corresponding boundary conditions are 
= 5.9 177.0 
T,=0 dS; f S,9~ d TS tale tre 479 dT! 9 (3.8) 


de tig ae art” fae 
aie (Seale 
This eigen-value problem was solved numerically to give y=71=y,, say, which is about 
twice y.=42. If the present problem is to be converging, y must tend to a finite limit on 
increasing the number of harmonics concerned. From this point of view, the above result 
is rather unsatisfactory. 

When harmonics up to the sixth degree are included in the argument and 2 or @ is made 
very large, we get 


S(6)=S;(E)= T(E) =0 (3.9) 
in the similar way as in (3.3)-(3.5). We have eight equations of the types as in (3.2) among 
eight non-vanishing radial functions 

TE)=2T(E), To(E)=xTe(E), T(E), S3°(E), THE), SHE), SPE), THE). VGAe 
Putting V,°=1—& and V.?*=(1—&)? as in (1.7) or in (3.7), we can reduce these equations 


into 


iT? 28 d’T,®  180d2T,° 183636 ; dT,° < 
ae Fae) tee ia xiaag ES —©) “ae +(45—658)T.* |= 


(3.11) 
The corresponding boundary conditions are 
5 dose 31 d°T,° aT, dT,° 
T’=0 , Sst + 26——*— + 154 = = 
6 dé + dEs a5 de? i dE 0 (1:2) 


at €=1. This eigen value problem was solved numerically to give y=117=y, say, which 
again differs much from yz and y, obtained before. 


Proceeding in this way, we can get the following general results. These results were 
obtained by Takeucut and Butuarp and reported partly in BuLLarp and GELLMAN’s paper 
(see p. 266-267) already cited. In the m-th approximation, say, in which harmonics up to the 
n-th degree are taken into account and 2 or @ is made very large, we get iS,") Se Ser 


20, soe, S84; Ti, Te,->:, Tn’, Te, Tes) >, TS =0 and we have eisht differential 
Bi accos » the types as in (3.2) among ite non-vanishing radial functions 


72, @)=2Th, ©» .Ta(6)=eTnl)y 


ie s(€); Dae s(€), fh ol o(&), eh i(€), per 1 (€), Tn?5(E). (G13) 
Putting V;°=1—€ and V.?*=(1—&)? as in (1.7) or in (8.7), we can reduce these equations 
into 
d*T,®  4(n-+1) dT," — 2n(2n+3) d2T,° 162n(n—1) 
+ roost = y 
d&* E dé* - dé (2n—3) (2n—1)?(2n+1)*(22+ 3) ed—é)" 


x | A(mE(1—€) a +| Bin)— we : 7.» |=0 ’ 
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A(n)= 4n'+12n3— 13n?—33n-+18, B(n)=8n' + 36n' + 1023 —1052?—632-+.54, 
C() = 32n° + 16490! + 10023 —485?— 41772-4306. (3.14) 
The corresponding boundary conditions are 
rT dSp_ ds Tao ha 
08 : + 2 —] Sree ® ee dT," = 
dé (2n—1)Sn_1~ dé 42(2n41)" op + 2(ma-+1) (2n—1)- dé =(0 ~ (3.15) 


at €=1. The eigen value yn, say, in this 
general case is shown to be between O((2m)°/?) 
and O((272)°/*) as moo, In the (72+1)-th approxi- 
mation, in which harmonics up to the (2+1)-th 
degree are taken into oe and _« or @ is 
made very large, S°_,(&), S%_,(&), T_,(&) and 
T?S (€) which did not vanish in the #-th 
“approximation become zero anew. Thus 
it seems that as m tends to infinity, every 
radial function vanishes and we have also a 
strange result that the eigen value y=4zKaV)2.2¢ 
increases without limit with increasing 2. 

The most fundamental assumption in carry- 
ing out the analysis in the present section is 
that in dealing with (3.13), for example, we 
are making the -th approximation for the 
one and the same system, the first approxt- 
mation of which is the S,°-T,°- T,2°- T.2* dynamo 
in section 2. From this point of view, yn 
obtained by solving (3.14) and (3.15) is the 
n-th approximation for the ‘‘true”’ y and yp, 
must tend to this finite limit as 72>, if the 
present problem is converging. The result in 
this section is against this anticipation. 

§4. Looking from a little different point of 
view, we can get very hopeful interpretation 
~ of the results in the last section. In dealing 

with (3.13), for example, we considered in the 

last section that we are making the m-th 

approximation for the system in section 2. 

From the present point of view, this inter- 
j pretation was wrong. Our new interpretation 

is in dealing with (3.13), we are studying the 

n-th eigen function (or rather eigen group of 
functions) corresponding to the m-th eigen 
value yn, which, from the present point of 
view, is not the m-th approximation for the 
first eigen value as was supposed in the last 
section. There are several reasons in _believ- 
ing that our present point of view is all right. 

At first, yn increases without limit as 2->09. 

This should be so, if yx is the m-th eigen 

value in the present problem. Secondly, only 


the (”—4), (#—3), --+, m-th harmonics belong 
to the group of fields corresponding to yn. 
This also should be so, if yn is the eigen 
value corresponding to the 7-th ‘‘ over tone’’. 
We are, so to say, in the similar position as 
those who are studying the m-th over tone of 
vibrating string streched between two fixed 


points. The differential equation for the 
vibrating string is 
ds 
yo=0, : 
dé? aa (4.1) 


where y is the eigen value to be determined 
by the following boundary conditions 

SE)=0 at €=0 and Ir (4.2) 
(4.1) is of the similar type as the equations 
in (2.1), say. If we try to get the m-th over 
tone of the string in which harmonics up to 
the 2-th degree are included bat not the (m-+-1)- 
th harmonic, we shall have 


y 2 Fi 
=("7) =Yn, say and S=sin ae =Sn,Say. 
(4.3) 


In this case, yn=O(*) as n.tends to infinity, 
to be compared with Yn = O((20)8/? )~ O((202)*/?) 


obtained immediately after (3.15). The m-th 
eigen function in this case is simply S= 


sin(z7z&/1), to be compared with the eigen group 
of functions in (3.13). This difference in 
simplicity is not so serious. By solving the 
eigen value problem (3.14) with (3.15), we can 
determine the relative magnitude of each 
radial function in (3.13). Calculating the field 
corresponding to each radial function and 
summing up the results thus obtained, we 
shall have the expression H(r7, 0,6; yn) or 
more simply H(yn), say, for the field corres- 
ponding to yn. This expression Hn) 1s, at 
least in principle, as simple as S, in (4.3). 
The physical reason why we get the larger 
eigen value y for the higher over tone is as 


follows. As is seen from (3.13), fields of the 
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higher harmonics, which are known to 
dacay more rapidly are contained in the higher 
over tone. Thus, as is shown in (2.4), in 
order to maintain the higher over tone, we 
must supply it with the more energy, and we 
have the larger y. From the above consid- 
eration, we can get rough idea on what the 
field H(yn) will become when it is driven by 
fluid motion for which y=4zKaV2,., is larger 
or smaller than y,. In the former case, given 
more energy than dissipated, the field H(yn) 
will develop, while in the latter case, it will 
decay. Suppose here that we have in the 
earth’s core y just equal to #2. Under this 
condition, the fields corresponding to H(y,), 
H(ys),-** will decay soon, even if they existed 
at a certain time while the field H(y2) will 
be sustained. This is very similar to what 
is occuring in the earth. Thus in our earth, 
magnetic fields of relatively lower harmonics 
(the dipole field, for example) remain almost 
constant through the geological time, while 
magnetic fields of relatively higher harmonics 
are changing quite rapidly. From the present 
point of view, we have no more thing to do 
for the self-sustaining dynamo model in section 
2. It is complete by itself. It is the dynamo 
corresponding to the first energy level, so to 
say, ¥, not the first approximation for some 
dynamo, as was supposed in section 2 and 3. 
The eigen group of functions for y, is shown 
in Table 1. It will be interesting to work out 
the eigen group of functions for higher y»™. 
A word must be added here on that we have 
T(E), for example, in both H(y.) and H(y,). 
There is no interaction through this (or rather 


these) 7,°(£) between Ay.) and A(y,). For 
example, in the above supposed case, in which 


we have y just equal to y2, the T,°£) field in 
Hy) group will decay according to the decay 
law of the group, independently of the T,°(£) 
field in H(y2), which remains constant through 
the time. 

Although the above results were obtained 
under the condition that 2 or a@ is very large, 
this condition is not so serious. This is because, 
except for small a, the eigen value y and the 
corresponding eigen group of functions vary 


very slowly with a. In fact, numerical results 
obtained by BuLLarp and GELLMAN (see Pp. 246 
of their paper already cited) show that there 
is almost no difference between the results 
obtained when a=(z/y) is 10 and . Thus 
we may expect the similar results as above 
when 2 or @ is not so large. 

It is now suspected that the earth magnetic 
field was reversed from time to time through 
the geological time. At present, we do not — 
know exactly how this reversion took place. 
From the results obtained above, however, 
we may construct the following model for the | 
reversion. Suppose that at a certain time in~ 
the earth’s history, convective fluid motion in 
the earth’s core became weak in strength and 
the corresponding y=4zKaV2,. became small- 
er than its stationary value y.. Suppose 
further that, after a time of smaller y, the 
convective motion and the corresponding y 
came back to what they were before. Through 
the time of this change, the rotational motion 
of the fluid was kept stronger than the con- 
vective motion, as was assumed in the present 
paper. This is an event not so unlikely to 
occur in the earth core. It is not necessary 
that any motion changes its direction. The 
convective fluid motion has only to make a 
certain fluctuation around its stationary value. 
This kind of thing, more drastic than supposed 
here, was supposed to occur in the earth 
mantle in order to explain mountain building 
cycles (see, for example, D. T. Grices, 1939). , 
Now, let us see what will occur in our model 
earth. According to what was said before, 
the field A(y2) will decay so long as y is 
smaller than y:, and will vanish in a time 
comparable with its free decay time. Given 
afterward y equal to its stationary value Yo, 
the field will be restored. Whether the new 
field, however, is in the same direction as 
before or not will be determined only by 
chance. Thus, in some case, we shall have 
the reversion of the earth magnetic field. 


Appendix 
Recently, S. Karo and H. Taxeucut (Geo- 


* see appendix 


“different from zero. 
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physical Institute, Faculty of Science, Tokyo 
University) have made the numerical integra- 
tion of (3.7) under the boundary conditions 
(3.8) and obtained the eigen group of functions 
corresponding to y,=71. Ina different meaning 
and in a different way from here, the eigen 
group of functions has already been obtained 
by BuLuarp and Geutuman. As was shown in 
(3.5), S,:° in this group must be zero. In 
BuLLaRD and GELLMAN’s paper (1954), however, 
this function is not regarded as zero. In fact, 
in Table 6b in their paper which is concerned 
with the same problem as ours, we find S,° 
As was pointed out by 
them (see p 258 of their paper), this illustrates 
a general difficulty in the numerical solution of 
complicated sets of differential equations by the 
finite difference method used by them. Anyway, 
this error may cause other troubles to make 
their result useless. From this point of view 
only, it is desirable to make the numerical 
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integration of (3.7) and get the correct eigen 
group of functions corresponding to y,=71. 

In order to compare our results with theirs, 
after finishing the numerical integration, we 
must transform our radial functions into 
theirs. The correspondence between their 
notations and ours is shown in Table 1. 


Table 1. 
BuLLARD and GELLMAN’S Ours 

7 & Ee 
Syd Yr) —&E” ey tune) 
Tr™7) — Ent Tn™E) 

A y 

E a=2/y 

Qs 6° Va (e) 
Q7/E eee Vile) 


We may disregard the difference of -+sign 


Table*2-a: 
r S29 T2/€ tise ® SEY aes Ty/é 
0 0.0 0.0 0.0 0.0 0.0 0.0 
"OL — 0.0%5000 0.071972 0.024971 0.045224 0.012662 0.06651 
0.2 — 0.032700 0.01580 0.03210 0.037283 0.0°5740 0.071123 
0.3 —0.01027 0.05341 0.08403 0.073113 0.071884 0.027205 
0.4 —0.06163 0.1291 0.1738 0.028834 — 0.01652 0.01160 
0.5 —0.1308 0.2671 0.3421 0.02221 —0.1052 — 0.05223 
0.6 —0.01842 0.5029 0.4456 0.04839 — 0.2136 — 0.3047 
0.7 0.3619 0.8115 0.022100 0.07627 —0.1141 — (0.7184 
0.8 0.6601 1.0000 — 1.0169 0.07425 0.2117 — 0.9842 
0.9 0.6328 0.7779 — 1.4867 0.03269 0.3737 — 0.7793 
1.0 0.5007 0.0 0.0 0.0 0.0 0.0 
a a eee 
Table 2-b. 

Se Per ee 0.0 0.0 0.0 0.0 0.0 
0.1 0.000 0.001 0.002 0.000 0.000 0.000 0.000 
0.2 0.000 0.005 0.011 (0.000 (0.000 0.001 0.003 
0.3 —().006 0.019 (0.034 0.001 ().000 (0.004 ().008 
0.4 —0.038 ().056 0.108 ().004 —().018 ().009 (0.020 
0.5 —0.089 0.165 0.287 0.013 — (0.087 —().017 0.036 
0.6 — 0.030 0.431 0.402 0.037 — (0.166 --(.157 0.020 
0.7 9.231 0.810 — 0.038 0.066 —(0.100 —0.418 — 0.068 
0.8 0.435 1.000 — 0.930 0.064 0.107 — (0.590 UC 
0.9 0.392 0.729 — 1.226 0.027 0.204 —().458 —().166 
1.0 0.236 0.0 0.0 0.0 0.0 0.0 —0.138 
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Table 3-a. 
Vi Sy; T2/é Die 
0.0 0.0 0.0 0.0 
0.1 0.06599 — (0.02287 — 0.05648 
OZ 0.04219 — (0.01506 —().02829 
0.3 0.1742 — 0.03324 —(0.03227 
0.4 0.4156 —().01179 0.06775 
0.5 0.4819 0.1382 0.2579 
0.6 — 0.03791 0.4625 0.2394 
0.7 — 0.9820 0.8374 — 0.2908 
0.8 — 1.6279 1.0000 — 1.0470 
0.9 — 1.6883 0.7394 — 1.2087 
1.0 — 1.5322 0.0 0.0 
Table 3-b. 

if 8,9 T2/€ tines 
00 00 00 0 

0.1 0.019 — 0.002 —0.006 
0.2 0.084 — 0.014 —0.021 
OS 0)239 — 0.022 0.001 
0.4 0.444 0.022 0.121 
0.5 0.388 0.199 0.273 
0.6 —(0.257 05382 0.169 
(OY —1.211 0.882 — 0.394 
9.8 — 1.764 1.000 —1.085 
0.9 — 1.741 0.717 —1.177 
1.0 — 1,568 0.0 0.0 


between their radial functions corresponding 


to the magnetic fields and ours, since in the 
eigen value problem as ours, we are not 
interested in the absolute values of the eigen 
functions. 

On the other hand, the difference of sign 
in the radial function corresponding to the 
velocity fields cause some difficulties. Any- 
way in Table 2-a, are shown the values of 
Butiarp and GELLMAN’s radial functions trans- 
formed from our calculated results. The 


corresponding results obtained by BuLLarp and 
GELLMAN are shown in Table 2-b. Comparing 
these tables, we see no serious difference 
between them other than in S,°. This is a 
very fortunate result. In order to get the 
eigen functions corresponding to 4, we have 
only to solve a differential equation of the 
fourth order. Compared with the rather compli- 
cated method by BuLLarp and GELLMAN to 
get them, this is a considerable simplification. 

By the way, we shall make the similar 
comparison of the radial functions correspond- 
ing to y,=42. In Table 3-a, are shown the | 
values of ButLarp and GELLMAN’s radial func-~ 
tions transformed from Table 1 in section 2. 


The corresponding values by BuLLarp and 


GELLMAN are shown in Table 3-b. Here also, 
the agreement between the results is good 
enough. 
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Abstract 


Wave generations from a line source of SH type in a superficial layer resting on a 


lower layer are studied. The movement at a 


point is shown to be composed of pulses, 


each representing a disturbance arriving at the theoretical time for one of the reflected 
and refracted waves. When we superpose these component waves, we get, at points far 


away from the wave origin, the theoretical 


movement which is very similar to that 


predicted by the normal mode solution. The reflected and refracted waves are shown to 


appear at the respective critical distances. 


Thus we can make clear how the transition from the aperiodic motion near the 
origin to the complicated oscillatory motion far away from the origin takes place. 
The most significant result obtained is that the amplitude of the reflected wave is very 


large at the respective critical distance given 
Examples which have some connections w 


§1. In connection with the generation of 
Love waves and the propagation of sound 
waves in the ocean, many theoretical studies 
have been made on the wave generations from 
line or point sources in a superficial layer rest- 
ing on a semi-infinite lower layer (H. Jerrreys, 

“M. Musxarts, K. Sezawa, S. Sato, C.L. PEKERIS, 
M. Ewina, H. Honpa). Formal solutions of these 
problems are obtained in the forms of com- 
plex integrals. The most successful approxi- 
mate way in evaluating the integrals is to use 
the steepest descent method or the stationary 
phase method. When applied to the present 
problem, it gives the dispersion formula in 
terms of the group velocity. Using the 
dispersion formula thus obtained, we can 
successfully interpret the wave phenomena at 
points far away from the wave origin. Being 
essentially the principle of complete wave 
interference, the steepest descent or stationary 
phase method is not useful at points near the 
wave origin, nor immediately after the time 
of beginning of the movement. At such 


lop (Mee 
ith this result are given in the last section. 


points and at such a time, another method 
must be used which is complementary to the 
steepest descent or stationary phase method. 

As such a complementary method, we shall 
in the analysis to follow use a kind of ray 
theory. We shall expand the integral in a 
series of exponentials, each representing a 
disturbance arriving at the theoretical time 
for one of the reflected and refracted waves. 
This method has already been used by some of 
the authors above refered to. But it seems to 
the present authors that its importance has 
never been fully appreciated. To make clear the 
importance of the ray theory in the problem 
of wave generations will be one of the objects 
of the present paper. In the case of explo- 
sive wave origin, the movement near the 
origin is very simple. On the other hand, as 
is shown by the steepest descent method, the 
movement far away from the origin is very 
complicated. Its oscillatory character seems 
to be due to the multiple reflections in the 
superficial layer. But it has never been 


22 Hitoshi TAKEUCHI and Naota KOBAYASHI. 


shown how this transition from the aperiodic 
motion to the periodic one takes place. In 
the analysis to follow, by means of the ray 
theory, this transition will clearly be shown. 

§2. As the problem of sound wave propa- 
gation in a superficial fluid layer is shown 
to be mathematically quite similar to that of 
SH wave in superficial elastic layer, we shall 
confine ourselves to the study of the latter 
problem. Refering to the rectangular coordi- 
nates in Fig. 1, we shall denote the density 


r ] 


ae ie yee pe 
fe) Rk’ 
v 
| ee Eaves 33] 
Fig. 1 


and rigidity in the superficial and the lower 
layer by 0, and 0’, yw’, respectively. The 
thickness of the superficial layer is H, and 
the wave origin is at (@#=y=0,z=D). The 
only non-vanishing displacement component 
of the SH wave is the y component. Denoting 
this by v and assuming 


poc ei(Mt +Ex) ‘ 


(219) 
we have the following differential equation 
for the z part of v, v(z), say, in the superficial 


layer 
a 
zn 
where 


B=2—R , 


™ B* )ov2)=0 (2.2) 


Bee 055 Pl ow ig ay 
et Vs 

V, being the distortional wave velocity in the 

upper layer. Similarly, in the lower layer, 

we have 


(42-8? )o@)=0 


pr Bk, 


Cl pee Pe 
k =)", Ve > 


Vs’ being the distortional wave velocity in 
the lower layer. The line source at z=D is 
expressed by 


(2.4) 


y=” eivt (RR) 
21 


ey t (edt 
age a aE B 

for zsD, (2.5) 
where 

R?=(D—z)?+2? . (2.6) 
The z part of v in (2.5) is 

p= ton . (250 
Correspondining to this v, we have some 


secondary wv as the solutions of (2.2) and (2.4). 
Thus the total v in the upper and lower layer 
will be 


1 
y= Ae 8?+ Bes? + —— et B-» 
1 e Fe at 28 
toe” (DEI). 
v= Ce fou. 2]. (2.8) 


where A, B and C are some constants to be 
determined by boundary conditions. The 
boundary conditions to be satisfied are the 
continuity of stress and displacement at the 
boundary z=H and the stress vanishiment at 
z=0. These boundary conditions are satisfied 
by making v(z) and 4“(dv/dz) continuous at 
z=H and dv/dz=0 at z=0. Having thus 
determined A, B and C, we get v(z=0) as 
follows 


; +o p-BD1 Ki p8(D-2H) be 
v(e=o)=oir| See ede, (2.9) 
where 

uB—W BE ; 
— pB+HB ae 


Putting D=0 for a while, and expanding the 
integrand in (2.9), we have 


EY 


Wave Generations in a Superficial Layer Resting on a Semi-infinite Lower Layer. 23 


co 
a! 


BL DV= SS Oe ¢ 
n=0 
n= ner | ri G EA ANMB+iEx dé : (2.11) 
where 


In order not to make v(z) infinite at zoo and 
in order to satisfy the radiation condition 
(SOMMERFELD, -A.: 1949, p. 189) at x>o, we 
must have 

the real part of (6) xghavel jel ss) (2.13) 
We can satisfy these conditions by interpreting 

the integral in (2.11), as the limit of the con- 
~ tour integral 


m= Ererrt| ee i 
L 
a A 4102) ye 
VCR 

when semi-circles around four points +k and 
+k’ are made vanishingly small and VW ¢?—R? 
and /¢2—,”? at the left end of Z are assumed 
to be positive real. The integral in (2.14) is 
of the similar type as that we had in a 
previous paper (TaKkerucut, H. and Kopayasut, 
N.: 1955). Thus following the way in that 
paper, we shall transform (2.14) into 


(2.14) 


EN kegs e 


Aes “aad EE elas sh oa 
a n en UV ORTH VY CR? 
ey ES ee) ge (2.15) 
V C2—R? 
where E->A-—D is the path on which 
Z=i€x—2nHV Ck? 


(2.16) 


in (2.15) is pure imaginary. The points EF 
and D will be made far away from A in the 
following analysis. 

§3. In order to transform (2.15) into a 
more convenient form for numerical calcula- 
tions, we shall make some mathematical 
considerations in the present section. 

(a) At first, we shall assume that the wave 
velocity in the lower layer is larger than that 
in the superficial layer. This means k’<k as 
is shown in Fig. 2. Now it can easily be 
shown that the coordinate € at A and Z in 
(2.16) there are 


So " , Z=—tkr, 1r=Vx'+(2nH) , 
(ail) 
respectively. 


(b) €, Ve2—Re and /¢2—R on AD are of 
the forms —(+)+2(+), where (+)’s are some 
positive real numbers. € on AF is the complex 
conjugate of that on AD, whereas ¢?—R? 
and Y €2—k? on AE are —(complex conjugates) 
of those on AD, respectively. 

(Cc) Z=1€x—2nHy @—R2 on AD and AE 
can be put 


Z=—ikrcosh¢, O0<@<+o, (38.2) 
where 
c= * (—x cosh ¢+i2nH sinh 4) , 
Z 
d€=+V C?—k do , Gd) 


respectively. 
(d) By means of (3.3), (2.15) is transformed 
into 


) 


m= 26 | V(o) exp {i(pt— kr cosh ¢)}d¢ , (3.4) 


where V is the real part of 
esta saga 
UV O—B+WV COR? 
on AD. 
(e) Putting 


1 


Cah =pbt, b= a, (3.5) 


and assuming S(t) instead of e’” in (3.4) for 
the time variations of the displacement and 
stress near the origin, we have 
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n= 26a) “Vie )}St—br cosh 66 86) 
0 
(f) The stress near the origin is given by 


Tees e0. Tae io — BS). (3.7) 


In the following numerical calculations, we 


shall assume 
S@=0 for <0, 


St=0)=0 , St=0)dt=finite, S, say. (3.8) 
(g) For the time variation in (3.8), we 
have 
tea ieee 
ag br sinh wy) cosh =t/br=V gt /r 
0<¢<~, (3.9) 
where 
C= : (—x cosh $+22nH sinh ¢) , 
x 
r= V 22+ (2nH)? (3.10) 


auetal WED) Te 
the real part of 


pate Ase ean on AD , 
UV E?7L EMV C2 — Mm 
(3.11) 
m being m=b’/b=V,/V,’. It is to be noted 
here that the real and imaginary part of 
Ve? or VE"2—m? on AD are negative and 
positive, respectively. (3.9)-(3.11) are the 
formulas by which we can calculate the dis- 
placement on the free surface for the wave 
origin (3.7) and (3.8). 
(h) When 
a ee e, 2nHk’ 
ook or Fys 


2nH V; 

Ve ie 
(3.12) 

the point A comes to the left-hand side of 

€=—k and there occurs the similar difficulty 

as we had in the paper above refered to. ln 

this case, following the way in the above 

paper, we have the following additional v» 


x> 


Vy (additional) 


“Sa : 2EnSU Porte 
a? tls y- efi ant E j 
by’ wv ‘(b/b y—é 
bx 
Ses b’y ’ 
t 2nH b 1/ 
Oe 


whese U is 
the imaginary part of 


| ee se 


BY OBL WV CR 
Ve2—k and Ye2—R in (3.14) must be in- 
terpreted as z[(b/b’)?— 
being shown in (3.13). 
$4. The time of beginning of the elemen- 

tary wave in (3.9) is given by 


— 7 _Ve+(2nH}? 
Vs 


gy? and (2-1), £ 


(4.1) 


Fig. 3 


7 
x—2nH tan A) 


Fig. 4 


As is seen from Fig. 3, this is the theoretical 
time of arriving of the wave reflected » times 
at the lower boundary. The time of beginning 


of the elementary wave in (3.13) is given by 


fe Teo, «—2nH tan@ 
V.con@ Wate . 
sin 9 = Vs (4.2) 


V;! 
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=1.0 


Toh aioate esl 


(c) 


(d) 


Fig. 


| call aoa LL Sun 


0.9)" 


EG (g) 
oy 
2.0 
1.0b 
O- ae < — re 
mle 
(Gin) 


Vie 


29 


26 Hitoshi TAKEUCHI and Naota KOBAYASHI. 


As is shown in Fig. 4, this is the theoretical time 
of arriving of the wave refracted m times at 
the lower boundary. From Fig. 4, it can easily 
be understood that we must have the condition 
(3.12) for the occurrence of the m-ply refracted 
wave. Anyway, we see from the above 
analysis that the movement at a point is com- 
posed of the multiply reflected and refracted 


waves. In order to see the wave forms for 


the impulsive origin (3.8), putting tentatively 


are 
SV, 


p= 0: ’ Ve=O08Vea ’ (4.3) 


we calculated the movement at «=20A. 

The results of the calculation are shown in 
Fig. 5 and’6. «In Fig. 5, the) reflected? and 
refracted waves of the m-th order are shown 


separately, the time ¢ being refered to the 
arrival time of the direct wave. 


In the movement for each ”, that before 


the time of the respective infinte displacment 


is due to the refracted wave, whereas that 
after the time is due to the reflected wave. 
By the condition (3.12), the number of the 
refracted wave at «=20H is limited to 7. By 
this and one more reason to be given later, 
only the component waves up to the seventh 
order arésshtown in Fig. 5, In Fig. 6, sis 
shown the result of superposition of these 
component waves. The symbols Rj, Rs, ---, 
D, Ly, Ly,--- in this ‘figure show the arrival 


times of the refracted waves of the first, 
second,--- order, the direct wave and the re- 
flected waves of the first, second, --- order, 
respectively. From Fig. 5 and 6, we see 

(a) The apparent periods of the waves 
before the arrival time of the direct wave are 
rather long. They are shorter, the nearer to 
the direct wave. 

(b) The apparent periods of the waves after 
the arrival time of the direct wave are rather 
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short. They are shorter, the nearer to the 
direct wave. The results (a)-and (b) are very 
similar to those obtained by Pexerts (1948) for 
the sound wave propagation in the ocean. In 
this case, the waves before and after the 
direct wave are called the ground and water 
waves, respectively, and it was shown by 
PEKERIS using his group velocity curve that 
the apparent period becomes very short im- 
mediately after the arrival time of the direct 
water wave. 

(c) The refracted waves after the direct 
wave are masked by the reflected waves, 
» whereas those before the direct wave are clearly 
seen in Fig. 6. The reflected waves are well 
separated from one another. In fact, we may 
easily construct Fig. 6 from the curves for 
the component reflected waves in Fig. 5. In 
section 1, we said that the ray theory we are 
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now using is complementary to the steepest 
descent method and is not suited to the study 
of wave phenomena far away from the wave 
origin. 

This remark now requires a little correc- 
tion. From what was obtained above, we 
see that our present method is fairly suited 
to the study of wave phenomena away from 
the wave origin. 

We shall now give the reason why we 
omitted the contributions from the reflected 
waves of the orders larger than 7 in Fig. 5 
and 6. As is seen from (3.9), the amplitude of 
the reflcted wave v, becomes infinite for ¢=0 
or at t=br (the arrival time of the reflected 
wave in question). This becomes so by the 
factor sinh @ in the denominator of vn» which 
is zero for ¢=0 or at t=br. Thus the value 
of vn for a very small ¢, dd, say, will be 
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br cosh(¢=0)do ; 
e’ in V(C’) in (4.4) may be taken to be equal 
to that for ¢=0, that is, €’=—(a/7) by (3.10). 
The correspondence between d¢ and dt, which 
is the time deviation fom ¢=07, is given by 


(4.4) 


t+dt=br | cosh(d=0)+ sinhg=0) doh 


400s (ad) 
1 fir 


ae Voor (4.5) 


or 


Thus we have the value of v, at t=br+dt as 
follows 


= V26nS _ Vie =—(a/n)} 
Vn(br+ dt)= V bHdt [(a/H)?+(2n)?] 1/4 7 
r=V 4 OnE. (4.6) 


Keeping df in (4.6) as a small quantity in- 
dependent of ” and 7 and omitting the com- 
mon factor S//YbHdt, we calculated the 
relative amplitudes of the reflected waves for 
several 2 and 2/H. The results of the calcula- 
tion are shown in Fig. 7. From this figure, 
we see that the amplitude of the reflected 
wave of the #-th order is very small up to 
the critical distance x=2nHV3/VV.2—YV;3?, 
which is, by (3.12), also the critical distance 
for the appearence of the refracted wave of 
of the m-th order. The existence of the com- 
mon infinite factor will make the above dis- 
cussion more complicated. But from what 
was obtained above, we may safely conclude 
that the number of the reflected waves at a 
point is the same as that of the refracted 
waves there. Thus we have seven reflected 
waves at «=20H. This kind of relation be- 
tween the refracted and reflected waves is 
familiar in optics. In optics, we have the 
totally reflected wave when there is no wave 
to go to the lower layer. In every text book 
on elementary physics, we can find illustrations 
of this phenomena for the ray going, for 
example, from water to air. Almost the same 
thing is occurrring in our present promblem. 

§5. As is shown in (3.8), the movement at 
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the origin is very simple, whereas as was 
shown by Pexeris, that far away from the 
origin is very complicated and oscillatory in 
character. We can now understand how this 
transition from the aperiodic motion to the 
periodic one takes place. It is due to the 
appearences of new waves by the multiple 
refraction and reflection. As we go further 
from the origin, we have more refracted 
and reflected waves and thus more complicated 
movement. It is our opinion that this point 
of view is very useful in understanding the 
wave phenomena in a layered medium and by 
using it, we can get from seismograms, for 
example, more informations on the under- 
ground structure. 

We shall now give examples which have 
some connections with the above thoretical 
results. Hughes (1949) developed a method 
to meassure the wave velocities through 
elastic samples. He used the supersonic 
wave and his samples were in the forms of 
bars. An impulse was given to the one end 
of each bar and the resulting movement 
was recorded at the other end. He obtained 
the more complicated movements for the 
longer bars. According to him, the movement 
is composed of the waves, the arrival 
times of which are given by 
Va'—Vety* 


L ( 
{=m nD~ 
” V> ar VoVa 


n=(0, 1, 2,---, m=1,-3, 5,+°-, 7ho.) 


where D, ZL, Vp and Ver are the diameter and 
length of the bar and the compressional and 
distortional wave velocities through it, re- 
spectively. #2 and » in (5.1) are subject to the 
restriction 


mL=nD tan 6 . (x2) 


The waves »=0; m=1, 3, 5,--- are the direct 
wave, the waves reflected once, twice, --- at 
each end of the bar, respectively. The waves 
m=1; n=1, 2,---are the waves which traverse 
the bar with the velocity Vr once, twice, --- 
at an angle @ relative to the normal of the 
boundaries 


sin = (5.3) 


Vo 
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and then go with the velocity Vp» along the 
boundary to the end where they are recorded. 
The detailed analysis of this problem will be 
given in a forth-coming paper. However, it 
can easily be understood that all the waves 
in (5.1) other than »=0 area kind of refracted 
waves, and the condition (5.2) corresponds to 
our (3.12). Thus the appearence of the 
successive refracted waves seems to be the 
essence of this experiment. Observing all the 
arrival times of these refracted waves, we 
can determine Vp and Vx more accurately 
than when we observe only the arrival time 
of the direct wave. 

Using explosives, H. Tate and others (1953) 
have tried to determine the crustal structure 
in the United States. They say that almost 
all the “phases” in seismograms obtained 
are spurious ones caused by accidental wave 
interferences. Rejecting the spurious phases 
one by one, they arrive at the phase which 
seems to be the only “true” phase. Accord- 
ing to them, this phase corresponds to the 
totally reflected wave from the Mouorovicic 
discontinuity. This totally reflected wave 
has very large amplitude at the critical 
distance similar to that in (3.12) (of course, 
we must put 2=1 there) and is observed 
between 80 to 130Km. The thickness of the 
crust obtained by them is about 35Km. and 
the wave velocities above and below the dis- 
continuity are 6.1 and 8.1 Km/sec, respectively. 
_ By using these values, we can understand the 
reason why we have no reflected wave up to the 
skip distance 80 Km. The skip distance is the 
term in the study of radio wave propagations, 
where we have the similar phenomena. We 
may also give the reason why we have no 
refleted wave beyond 130Km. Thus, as is 
shown in Fig. 7, the amplitude of the first 
reflected wave decreases from the maximum 
at «=2HV:;/V V;’?—Vs? and becomes zero at 
about z=3.5H. This may correspond to the 
disappearence of the reflected wave at 130 Km. 
from the origin. In this connection, it will be 
interesting to study the existence of the second 
maximum in the first reflected wave and the 
appearences of the reflected waves of the 
higher orders predicted by the present theory. 


Anyway, by studying the waves totally re- 
flected at the Mounorovicic discontinuity, we 
are sure to know more about the crustal 
structure. We can also expect to know more 
about the underground structure by applying 
the similar method to exploration seismology. 

The hydrophone and_ sono-radio buoy 
method is reported by M. N. Hitt (1952) to be 
very useful in studying the crustal structure 
under the ocean. Typical records obtained by 
him are shown in Fig. 5 in his paper. They 
are composed of the waves which are separated 
from one another quite well. He says that 
the number of pulses for any one distance 
between shot point and receiver is limited. 
The number is the more at the more distant 
point. The path of each wave is shown in 
his Fig. 10. Each wave is called by him the 
multiple refraction wave, and some reasons 
are given why it is not called the multiple 
reflection wave. But it seems to the present 
authors to be rather difficult to separate the 
refracted wave from the reflected wave of the 
same order, the critical distances for the 
appearences of these waves being the same. 
As is shown in Fig. 5, each refracted wave 
comes earlier than the corresponding reflected 
wave, but its amplitude is so small that it is 
liable to be escaped from our notice. In fact, 
as was said at (c) in section 4 and was shown 
in Fig. 6, the refracted waves after the 
direct wave are masked by the reflected 
waves. From this reason, it seems to the 
present authors to be better to call the waves 
in Hiw’s Fig. 10 the multiple reflection waves. 
They are the reflected waves from the lower 
layer of varying velocity and analogous to the 
reflected radio waves from the ionospheres. A 
relation between the total reflection and the 
dispersion formula has been pointed out by 
several authors (e.g. Toustoy, I. and Uspin, 
E. 1953). The appearence of totally reflected 
waves of the higher orders at the respective 
critical distances to form the typical oscillatory 
movements at the point far away from the 
origin is probably another aspect of the 
above relation. 

The above analysis has been made for the 
SH wave. The corresponding problems for 
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the P-type origin will be made in a forth- 
coming paper. However, we are sure that 
the essence of this kind of problem is exhausted 
in the present paper. 
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On Microseisms Produced on the Surface of the Earth due to 
Passage of Storm over the Deep Sea. 


By 
Sushil CHanpra Das Gupta 
Krishnagar College, Krishnagar, West Bengal, India. 


Summary 


Several attempts have been made to explain the microseisms, the latest being that of 
Longuet-HIGGINS based on the idea that they are produced by stationary waves. In this 
paper the author has attempted to show that in case of storms, microseisms are also likely 
to be generated by progressive waves. The character of microseisms has been proved to 
be of the RAYLEIGH type approximately. The average velocity and wave length of micro- 
seisms have been shown to agree approximately with the observational results referred to 
in this paper. The period-relation observed by BERNARD and DEACON viz. that the period 
of microseisms is half the period of sea-waves supposed to generate them, is also shown 
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to agree. 


has been shown to depend on the geological structure and not on the distance. 


The ratio of horizontal and vertical maximum amplitudes has been found and 


In conclu- 


sion the possible reason for absence of microseisms on the approach of monsoons has been 


given as noted by GHERZI. 


Introduction 


Waves produced on the earth due to distur- 
bances at sea have been considered by several 
authors. GimMore (1946) and Ramirez (1940) 
particularly made elaborate observations 
regarding the seismic-waves due to passage of 
storms over the sea with comparatively high 
velocity. Others, such as GuTENBERG (1931), 
BAnERgI (1930) and Guerzi (1923; 1928) obser- 
ved what are called micro-seisms due to vari- 
ous types of disturbances on the sea. BANERJI 
found microseisms due to the coming of mon- 
soon in the Arabian sea on the Indian coast. 
- GHERzI observed on the China Coast micro- 
seisms due to the passage of typhoons but 
none due to monsoon. Among the many 
observations Ramirez observed the average 
velocity of microseisms to be 2.66km_ per 
second and the average wave-length as 14.76 
kms. Average period of microseisms was 
observed by him to be 5.2seconds and the 
average ratio of vertical to horizontal com- 
ponent of displacement to be 1.21 at St. Louis. 
GitmorE observed that amplitudes increased 
from 17mm to 30mm when the wind velocity 
increased from 55 knots to 85 knots per: hour. 


Microseismic activity at Richmond decreased 
from 32mm to 20mm as the strom passed 
over Florida when the intensity of storm 
dropped from 85knots to 40 knots. When 
the storm reached the land, microseismic acti- 
vity decreased very much. Observed periods 
of ordinary sea waves caused by hurricane winds 
average about 10seconds with wave-lengths 
of 500 to 700 feet. Velocity of sea waves sel- 
dom exceeds 30 miles per hour. 

BeRNARD (1946) obtained evidence which led 
him to believe that microseismic oscillations 
have periods which are half those of the sea 
waves which give rise to them. The same 
ratio was observed by Deacon (1947) between 
microseisms recorded at Kew and sea-waves 
recorded at Perranporth on the north-coast 
of Cornwall. 

About the origin of the earth-waves of the 
type of microseisms, there were three theories 
in 1937. The first was that of GUTENBERG 
who considered that these microseisms origi- 
nated in surf on rock-bound coasts. The 
second was that held by GuerzI, BANERsI and 
others that they were generated at the centre 
of storms at sea. The third opinion was that 
of Braprorp, Ler and others that they had 
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their origin in certain types of frontal distur- 
bances anywhere on land or sea. Discrimi- 
nation between these possible origins by fur- 
ther statistical correlations was hopeless, 
because all three often existed simultaneously. 
Recently Longuet-Hiceins (1950) has given a 
theory of microseisms considering the MIcuHE 
(1944) term of pressure-relation in a fluid. He 
has found that the total transmission of energy 
to the sea-floor is of the right order of mag- 
nitude to account for microseisms. He con- 
sidered that stationary waves give rise to 
microseisms which have period-relation as 
observed by Berrnarp. DarpysHIRE (1950) 
considered three storms to show the possibility 
of stationary waves being generated in an 
attempt to substantiate the hypothesis of Lon- 
guet-Hiccins. But there is always a great 
possibility of progressive waves being genera- 
ted in the open sea due to storms and there 
is likelihood of generation of microseisms due 
to them. 

In this paper, the possibility of generation 
of microseisms due to progressive sea-waves 
produced by storms has been considered. 
First, it has been shown that the approximate 
period of microseismic vibrations is half the 
period of the sea-waves generating them. 
Secondly, the character of microseisms is 
shown to be approximately of RAYLEIGH type. 
Thirdly, it is demonstrated that the velocity 
of propagation is much smaller than that of 
RAYLEIGH waves and that the velocity depends 
mainly on the geological structure. Fourthly, 
the wave-length has been found to conform 
to observations. The ratio of the horizontal 
to the vertical amplitude has been obtained 
and lastly an explanation of the phenomenon 
that microseisms are observed in case of ty- 
phoons but not in case of monsoons has been 
given. It has been shown that compressibility 
of water is a reasonable cause for propagation 
of surface disturbance to the sea bottom. 

1. When the compressibility of sea-water 
is taken into account the velocity-potential 
function @ at a depth z is given, as shown 
by Wutpp.e (1935), by the relations 

@=exp (¢hct—tlx)[27 2 cosh (22) 


—(p?—277) sinh (uz)je-” alain) 


when Pore 

and ei ae 

and @=exp (tpct—ilx)|27v cos (v2) 
—(p*—2y?) sin (vz) Je7?” (1.2) 

when Poy 

in which yv=in, w=p—P—;? 

and a : 


We have the frequency equations in case of 
(1.1) and (1.2) as 


Sarees Ses 2 


LL 
—coth zh=— ‘ ik) 
r 0 Pp 2 oe 
yp P—y?—y? 
— cot vh= ; 1.4 
ey 2 i 
where hk is the depth of the ocean. rh is 


small, + being of the order of 2x10-§cm7} 
and-/z of the order of 10°cms for the deepest 
oceans. Equation (1.3) regarded as an equa- 
tion in 4 has one real positive root which 
approximates to /—y when /h is large. Hence 
neglecting the square of velocity and remem- 
bering that eet. (Ge) where 7 is the de- 
g ot z=0 

pression of the sea-surface from the mean- 
level, the pressure at the bed of the sea can 
be expressed approximately in the form 

P=goh-+-gon}cos Bip ect vn 

2rv 

or 


goh-+-gon cosh ie 
2ru 


~ sinh en} ant Kis 
We thus see that there is a finite disturbance 
of pressure at the bottom and is proportional 
LOM. 


Let n= Aei(ée—pet) | (1.6) 


where c is the velocity of sound in water and 


ae Hence the pressure disturbance at 


the ocean-bed is given by 


Cc 


Ago}cos vh— amt 168 sin vha|etee-vep 


2rv 
or 
Ago{cosh uh— Pats she }ote- vat : 
2r ue 
Writing be=p, 
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we have, ; Hence ¢,’ is a particular case of ¢, where 
Ago {cos ee a ae ae vaheteen %.=$1 and A,=0. Therefore, 
ary D,= Df + D,/" , 
or 
p2/e—27? D,= grad ¢» , ao.=0, 
1 ame . a 
Ago {cosh ah— ave sinh phen : O.D/’=0, div D/’=0 
aie, (cf. Havasu (1942)) . 
2. The equations of motion in the elastic 3. Let us take the x, y plane of the rec- 
earth are given by tangular co-ordinate axis on the undisturbed 
3 : P _ @D surface of the semi-infinite homogeneous iso- 
@ grad div D—b*rot rot D= or ’ (2.1) tropic elastic medium of the ocean-bed, z-axis 
rene being directed inward to the medium. We 
where a= ; assume that the disturbance travels with uni- 
form velocity v in a direction which is taken 
and pa | as a-axis. We are interested only in solutions 
p where x and ¢ appear as a combination of 
Taking divergence of (2.1) we get x—vt. A formal solution of (2.4) is given by 
lindiv D=0- (2.2) Er-etideriBy-tavt 
1 0? ee 
where a=A— a OF : where a=y/(1- reas ; (3.1) 
Again taking rotation of (2.1) 
: 3 and a, 8 may take any real value. Therefore, 
fla-rot D=0', (2.3) . 
2 = —qzt+iaxr+iBy—iavt 
where fa g , is {\" a, mes GEE 
b? Of 
aS Similarly the solution of (2.5) is given by 
Vet: D=grad ¢+rotA. é 
Therefore, A=\| ALG Be Veiner? alan 
nN) ’ res a : 
and (grad div—A)A=0. uae d= (1 )aee | Va 
‘Then = ites , Hence 
and A=A,+A, ’ feos) 
where Adi=0, Oab2=0 (2.4) p=\\"" Dia, Bje- ae site any sent fee 
and grad div A,—-AA=0, MEO. (yD ee a. 
D, : e-Vztiae+iBy—iavt ud B : 
Therm D=D,+D,, +(\" ket, #) 
where D,=grad ¢:+rot A, , hare Dwe=iad, Dy=i8d, De=—ae, 
= t As, 
' te eid eu ‘ (Diz=x-component of D,, etc.) 
div tls Be o ef et Dyz=tBAsztq Ary , 
rot gis ee a iV, ey A 
so that i= grad ¢, Diet pA 
and Ag =0. : : 
e2D It is clear that zaD,2+1%8D.y—qD2.=0. There- 
From (2.1) therefore ae (he fore, we may put 
ap? Doz=Q 1 , Dyy= q $2 
2 
so that ° “01 <0 


OF and Dz.= tad, +B») , 
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Let the boundary conditions be X,=Y,=0 and Z,=f(x, y,t) at z=0. Therefore 


Xz _ODz OD: _ 


ac Oe ade 0 at z=0 
or, [[ 7 —2iags a2 +a —asp,jermerr-idada=0 
Hence 2iagg+(q?+a?)pitaBd,=0 . (3.3) 
Next, we know MG tg 0 at 2=0 
So that [| - 28a —aben G+ Bg eee dads =0 
Then 218qo+aBdit(q?+B?)p.=0 . (3.4) 
Goneidering the relation an div D2 = at z2=0, we have 
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D > 
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where Da, B)=daa' (at +)—{(2— 7 Jar papel : 
Hence, w, v, w at z=0 are given by 
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v Se = aq tax+iBy—tavt 
Wor Nor Deena Be dadp . 


Since it is well known that uniform hydrostatic pressure has little influence on the seismic 
wave motion, we take the boundary conditions as 


—Z,= Ajet(fe-rit) = F (a, y, 2) , 


cos vh— Py oe a sin vh , 
v 
where A,=Ago * 
cosh nh— pille*— 27" -sinh xh 
aru 
from (1.7) when vi<ia<h+ut, Oa <1: 


so that the disturbance is proceeding in the simple form of a rectangle of sides 7 and hh. 
As two sides of the disturbed area are likely to lengthen in the direction of travelling of 
the disturbance, we take it to be approximately a rectangular one. Z, is taken to be zero 
outside this area. Therefore from (3.5) and (3.7) 


— g(a, B=a\\"" re, Y, Le *&(x-vt) -1BY 4 —vut)dy : 


Putting x—vt=y and f(x, y, t)=Aetee-71, 
tBL)(] —e-*Ca— £)11) 


rely ok. ae 
= = i(EN-—p t+Evt) p—ia@n—-ipy =A et(Ev-py)t ee a (3.8) 
g(a, B) |’ A etes Fi e dydy= Ai Bena 
Hence the displacements at z=0 are given by 
Aw? . +00 
(#)z-0= — Ta aa \\" 
TEN et 2 =f —i(@—£)11,)pia(a—-vt)+ipy 
a 1— 7) + Be iB1)(] —e-t@-Di)Je 
/( z dads , (3.9) 


j [ a(ar-+ 8 y/(a— ae 84) (1 ar4 p2— (2-5 "3 ,)a° soe" ‘le _a)p 


= _Au i(Ev—p,)t es 
(4)2-9= PPT ) 


—oo 


ol $)sae-a/ Bron |= Bprsr|acranaaeiemneiin 
a [Hees ary/ (1-2 ee ary/(1- ac +@— -{(2- Fes anpe0 Aa 


x dadB , (3.10) 


(v)z-0= — JM en NN 


2 2 2 a2 ( 1—e-#) 4 — et £)11)gtaCa—vt) +iBY 
AS at a?+28 —2y/ ae ‘(a +8 e/ (1 a+ r/ e- be 


[ees py (1.— aren /(1— Sarre {(2- 5 ‘ae.+ 26" fe —a)8 


xdadB. (3.11) 
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4. The expression within the sign of integration will be shown to be a function of the 
distance of the area of the storm from the position of observation. Hence, the period of 
disturbance of the waves travelling through the ocean-bed to the position of observation 


But as it has been shown by several observers that the 


2n 
will have the period 
EOP 
period of such waves is half the period of waves on sea due to storm, we get Ev—pi=2p1 


Of, V= oe 1 We therefore expect that the velocity of storms is three times the velocity of 


5 


waves generated on sea. From theoretical side, Jerrreys (1924; 1925) has shown that the 
velocity of wind will be at least three times the velocity of sea-waves. According to 
GinmMorE’s observations ordinary sea waves caused by hurricane winds average about 
10 seconds with wave length 500 to 700ft. The velocity of such sea waves is ordinarily 
30 miles per hour. Hence if we assume the velocity of wind as three times the velocity of 
generated sea waves we have, 


E 6°28 x 132 : 
=—=- X 2xryv=—__—=1° 
Ev On TV 500 16 approximately, 
Spi=3- ae -20 af AE a1'8 approximately. 


Thus an approximate agreement with observational result is found. In the absence of larger 
number of more or less exact data available, agreement could not be shown with greater 
precision. 

5. In order to integrate (3.9) we resolve the expression (1—e-‘®’)(1—e-"@-1) into 


1+ e- #8’ e-t(@- 91 — g-tBl__p-i(@-£) ly (5.1) 


which gives us (3.9) as consisting of four integrals. The first integral of them is 


Atco ae/(1—-2 da 4g ela (a-vi) +ipy 
i : 


[Heotonl (B ee 0 Be [eRe saef econ 


Substituting a=scos ¢, B=ssing, x—vt=r, cos 0,, y=n sin W;, the integral becomes, 


£1 [st 9 = > » 
\ | COS= 9VY1—7, cos? gy e'"18 cos (@-@) 


-2{4/1—r, cos? gY1—r; cos? y—(2—r, cos? y)?}(E—s cos g)s sin aoe ; 2) 


2 


Vv v 
where t,= = and. t= pe Now, 
oe b 


BO CORE ee gee 2 >t" Jn(718) Cos n(Y—B8) . (5.3) 


Substituting (5.3) in (5.2) and taking the first term i 
} of (6.3 i 
ik a (5.3) as a first approximation, we get 
i he cos? YY 1—7, cos? g Jo(7s) 
aves cos* YV/1—r, cos? 9 —(2—ry cos? g)? (Es cos Os Sine ; Oe 


In order to integrate (5.4) with respect to s, we choose a contour as in the following figure 
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with indentations at the origin and the point given by é -. Hence 
cos 


taking the principal value we get (5.3) as 


: cos? WIG cos? OD. en ) 
op eid ee ) 


~* {41/1 —r, cos? Soh Le cos? lites: cos? ¢)?} =. -sin g 


(5.5) 


the contributions due to the other pole being cancelled by similar 
expressions from each of the four integrals arising out of resolving (1—e-#!)(1—e-'@-4), 
In order to integrate (5.5) with respect to 9, we put 


cos e=5(2+7) and sin e=3,(2-—) where z2=e!?. (5.6) 

2 Zz 2% 2 
The contour of integration is the unit circle with necessary adjustments to be done from the 
following considerations of the poles and branch points of integrand of (5.5) so transformed. 


For the ocean bed where the ocean is deep, we take the mean value of the ratio Aches 
Vee 
as 3°5 at high pressure. We take then in (5.5) ies 3 = say. When the substitutions 
T. 


(5.6) are made in (5.5), the expression within the bracket in the denominator assumes the form 


cape g fap an 
Hence 


ae i 
=3(ee fsa arson eros DT] 


We are to find the zeros of (5.8) which lie ee dN 4 


42>=r —y 
within the circle of unit radius. We shall 2 T. 5.10) 
consider only the principal value of the inte- pal aes T6535 
gral (5.5). The expression within the bracket Seam darks Ele eee 


i . ] real roots lying within d 
os aa aetna pi The contour is then chosen to be the unit 


ee : 1 
TEESE SaaS a ota aPPFO-  Gircle with necessary modifications to exclude 
Saito Ves % : the C b 1 
ximately where umy/ 4x °95 corresponding to the singularities and branch points by Hen 
a method. It is clear that k>h>m. Taking 
X11, (5.9) then the principal value, the contribution due 


ati=+ 
: to the singularity +— to the integral (5.5) is 
v1 


na being small. Corresponding to sin¥ 
1 it . ° * To1" 
——) in the denominator of (5.5) there niy/1- 
=A 2) ast! (en), (5.11) 
are two singularities at +1. The branch mae EF me Nes =n) ay 
V1 V1 


points on the real axis are given by 


38 Sushil CHANDRA DAS GUPTA 


where F(z) is given by the expression (5.6). 
Similarly the contribution due to the singu- 


larity ee will be the same expression as 


C1 
in (5.11). Hence the total contribution will be 


ee eas, 
Fagin) 


Taking the next term of the expansion of 
(5.2) the contributions will be cancelled due to 


the two singularities oe . Hence there will 


V1 
be contributions only from the terms Jo, Js, --- 


of the expression (5.3). 

6. If we take similarly the term e-’®’ and 
write cos #,.=x—vt, 7.sind,=y—l, we get 
on integration of the second integral of (3.9), 
the same amplitude term as in (5.12). Similar 
are the other two integrals of (3.8) with 

%3 COS $3=x—vt—l, , 738in 03=y, 
1, COS O,=x—vt—1,, msin@s=y—l, 
amplitude term being the same as in (5.12) 


with the factors e’41 and Ga é nm), n= oF a) 


respectively. Hence Eoeibetion to w at a, 
from these two singularities is 
vay) 1- Tk v7 
tA? * NS eRe ae (7 
Ais 


ADE 


ye) 
-e fs) -1(% for 


approximately as x, is large. It is clear that 
the approximate velocity of propagation of 


the disturbance is SUP As it has been 
2E |x 
shown that €v—p,=2p, or v=3p, we have 
I,  v _1/4x 95 
the locity as ——~= VEX 40 > Be 65 
velocity as 2E lm, 3 . S00 bu 
A+2n 


approximately. For the medium where 
Lt 
=3.5 we take the value of the velocity of 


S-waves as 4.4km, per second. Hence the 
velocity of microseisms comes out to be 2.86 
km. per second. Ramirez observed the average 


velocity of microseisms as 2.67 km per second. 
There is thus an approximate agreement. 

The additional waves generated which are 
obtained from the singularities at +1 have 
low velocity and are not likely to be detected 
by observation. Hence the only waves to be 
recorded at a distant station are those given by 
(6.1). 

7. The wave-length is given by 


Die LE ox 
Sela EC ay ee 
For ordinary sea-waves caused by hurricane 
winds the wave length as observed by GILMoRE 
is 500 to 700 ft. when the period of sea-waves 
is 10seconds. Hence, (7.1) is 


D7 le OP ASSOR 2c gel 96 
~s- ? = —x-—xb. 
E Sig v/b2 €E SS een 
: 3p; 27 
Putting v= E and —-=10secs., we get the 


1 
average wave-length to be 3°3xb. If b=4.4 
km. per second the value of wave-length of 
the vibration is 14.6km. approximately, the 
observed average wave-length being 14.76km. 


8. For maximum amplitude, it is not a 
worthwhile attempt to find it without very 
exact data as it depends on the depth, velocity 
and the period of the waves. At least, it can 
be seen that the amplitude of the disturbance 
is of correct order. 

The nature of the roots of the equation 
(5.7) has been discussed in datails by Hayasur 
(1942). If observations are taken in this light, 
constitution of the earth below the ocean can 
be determined at different place. 


9. The horizontal component of displace- 
ment can be obtained similarly from the ex- 
pression (3.10) and (3.11) as in article 5. It 
is clear that, if we follow the method of 
article 5, the terms containing J) will be 
cancelled due to the two singularities corres- 


: uN ; : 
ponding to z+ : =-ta, given in article 5. 


Hence the first terms of the integral (LO 
and (3.11) will contain a(% 


Vy 
distance of the point of observation from the 
approximate extremities of the rectangular 


r) where 7 is the 
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disturbed area on the sea. Then the ratio of 


: 1 ; 
the rat h a ay 1eGedl. 
the maximum vertical to the maximum _hori- aa 3:5 3 1. This value 


zontal displacement is given by is dependent on v viz. on the geological struc- 
PRR Gee ture. Average value of the ratio was obser- 
Bb? ji- si ved by Ramirez to be 1.21:1 at St. Louis. 


2{—(2— 722) 4ay/1— HE ty) 1 Par} If the ocean bed is not homogeneous but has 
4 4 4 a slight variation of elastic properties in the 

(9.1) upper part, there may be slight dependence 
approximately. As 2, is the solution of (5.7) of ratio on the period, wave-length etc. But 


in terms of ch. we have the mean value will be mainly dependent on 
® geological structure. 
ay/1— ve ney/1— es = 3(2- peat) ; 10. The expressions for the displacements 


u, v up to first terms are given approximately, 
Substituting this i in in @. 1), we get the value ¢ of for 2, large, by 


ae 9 picer—vpt bos ve sore >) Su vee “he, in n\—1(% r:) 
1 


< D pacgv—vy)t 2 pe my ' ai aaa aC n)—3( r:) 
EL eae) 


In evaluating w, v, e“?-® has been substituted for cos(g—@) and real part of the integral 
has been taken. Hence the direction of horizontal displacement gives us approximately the 
‘direction of the disturbed area. The position of the disturbed area, therefore, can be deter- 
mined by means of two observing stations and observing the direction of horizontal displace- 
ment. The distance between the stations being known the position can be located. Taking 
more than one observation of the position of disturbed area the direction of motion of the 
storm can be determined quite accurately. 


; : 2& 
11. When the end points of the disturbed area are at very great distances viz. aa (n= 
2 2& 
1, 2,3, 4) is large, contribution to displacement due to rn = 1) or H( re ro) (n=1,273; 4) 
rs 1 o 
becomes too small to be detected at a large distance as 
2a 1 (4n?—12)(4n?—32)- - -(4n?—47—1°)| 
he~() cos (2-5 1s {13 (— mt ieze J 
: ie (4n?—1?)(47?—3?).- -(4n?—47—3?) 
+s1n (2G re (Op 1) f= 827-1 


when |z| is large and jargz|<x (Cf. Wuitraxer (1915)). So it is likely that the monsoon 
which approaches with a very large front is not likely to give appreciable microseisms. In 
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fact Guerzi observed that on the China-coast microseisms were observed when there was 
typhoon, but none with the coming of monsoon. 
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Nomographs for the Phase- and Group-Velocity of Love-Waves. 


By 


Yasuo Sat6 


Earthquake Research Institute, Tokyo University, Tokyo 


Summary 


Nomographs for the phase- and group-velocity of Love-waves were made. 


Such nomo- 


graphs were already made by the same author a few years ago, but the present ones are 
more convenient for the practical work of seismometry because of the range of para- 


meters and the precision obtainable. 


Nomographs convenient for the computation 
of phase- and group-velocity of Love-waves 
were made. Similar nomographs were pre- 
sented by the same author few years ago 
(SaT6, 1953). However, they were not appro- 
priate for the practical work, because a wide 
range of parameters was put into a limited 
space, and consequently the precision was not 
high enough for those engaged in observatio- 
nal seismology. 

In the present paper, in order to make the 
precision higher, the domain of parameters 
was restricted to those as may be found in 
the practical study of seismometry. 


. 
Sie = = (5) 


[4] 


Fo 


6(V/Vg) 


Fig. 1 


§1. Fundamental formula for the compu- 
tation of the phase velocity of Love-waves of 
the lowest mode is as follows (SatTé, 1951, 
1953): 


w =sec 9-[Arc tan{z sec 9V/ sin? 6—sin? 4}] 
Ae 


where 
w--+--frequency of Lovge-waves. 
xy =n / (wand x are the rigidity of 
the lower and upper medium re- 
spectively). 


cosec p= Vi/V8 (Va and V;* are the velo- 
city of S-waves in the lower and 
upper medium respectively). 
cosecO=V/Vs (V is the phase velocity of 
Love-waves). 
By modifying the above expression the 
following form is obtained: 
T (=period measured by the unit H/Vs) 


= 270 


ve a 
Vs 3 Vs Ve _- Ve 
(es v * Arc tan v7 % XV ya yral 
Ve2 V3 
Somat (ED) 


where 

HT is the thickness of the surface layer. 

§ 2. Explanation of the figure and direc- 
tion for use of the first nomograph. 

1” [1] is the axis of V./Vs. When acer: 
tain value of V;/Vs’ is given, we take a point 
A on the axis corresponding to this assigned 
value. 

2° wis the curve of 1/1/{(Vs?/V2)—(V#?/Vs’2)}. 
At first we choose some value of V/Vs. Follow- 
ing arrow line, we get Con [2] axis, which gives 
1/V (sin? 0—sin? 09) =1/V/{V2/V?)—(Vs/Vs")}. 

3° [3] is the axis of ’/#. We take a point 
D on this axis and then obtain the point & 
on the [4] axis. E gives 1/xV/(sin? 0—sin? 4). 

4° Finally the curve 8, which must have 
the same parameter V/V; as the curve a, 
gives G(T). (T=2z/o). 

Ranges of the parameters are: 

VPV "—0.65~0.97, 


V/s 
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V/V:—1.01~1.10 (every 0.01), 
1.10~1.46 (every 0.02), 
HE a Ba! Se 
T —0.8~9.0. 

§3. For the computation of the group 
velocity by means of the second nomograph, 
we give V;/Vs’, V/Vs and y’/u. The funda- 
mental formula is (according to SaTé, 1953) 


U det AD) 


Vo (V[Vs2+ Sx, Y)’ 


U,V---:: group and phase velocity of 
Love-waves respectively, 

2 Y=FOY) +14); 

F(Z)=1/Arc tan Z-(Z+1/Z), 

Y=Z/z=sec 0-1/ (sin? —sin? 9») . 


In the following figure (Fig. 2) the process 
of the calculation is illustrated. 


fo) {9 its)} 
E(V/Vs) 
U/V 
NI ee 
P(U/Vs) 
fe Al(Vs/Vs) 0 K (7) 
J 
§(w'7m) 


SATO 


1° As before we take a point A on [1] © 
axis which gives Vs/Vs’. 

2° + is the curve of 

1/Y =cos 6/1/ (sin? 6—sin? 4). 

If we choose a curve y with a parameter 
equal to the given value of V/Vs, then the 
point Hf on the [6] axis expresses 1/Y. 

3° Next, [7] is the axis of S(, Y). Since 
the parameter of the curve d in the second 
nomograph (right below) is y(=’/#), we can 
obtain the value of S(z, Y)=V {F(xY)-(1+1/Y*)} 
by combining 1/Y and x which is naturally 
determined when the distribution of matter is 
given. 

4° As the third step we obtain U/V. In 
the third nomograph (right upper), the curve 


¢ (V/Vs) gives 


U/V={1+Sx4, YIK(V/Vir+Sx(x, Y))- 

5° If we want to make clear the ratio of 
the group and the phase velocity of Love- 
waves, the computation is finished by the 
above procedure. But if we want to make 
clear U/Vs, we must multiply the value thus 
obtained by V/Vs by means of the next 
nomograph (left upper) and get the final 
result. 

§ 4. The precision of the results obtained 
will be high enough for the practical work. 
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Chemical Equilibrium within Self-Gravitating Planets 
and Internal Constitution of the Earth 


By 
Yasuo SHIMAZU 
Institute of Earth Sciences, Nagoya University, Nagoya. 


Summary 


Chemical structure of a simplified model of the planet of which the components are 
Fe-Si-O is discussed. Chemical equilibrium and hydrostatic equilibrium are assumed and 
emphasis is placed on gravitational separation due to density differences and chemical 
interaction among components, which control the equilibrium arrangement of stable phases. 
It is concluded that the present gravitational field of the earth is too strong for reasona- 
ble arrangement of the components. 

Models with probable gravitational field are discussed and the following conclusions 
are obtained; (a) in 2,000°K isothermal model, Fe.SiO, content shows a maximum at a 
certain depth and SiO, and FeO which are at shallower and deeper depths respectively will 
disappear soon; (b) in 2,000°K (surface) 6,000°K (center) model, the radius of the Fe- 
core is one-half of the earth’s radius. The mantle is composed of Fe.SiO,, SiO, and 
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isolated O. 


1. Introduction 


In two previous papers, hereafter referred 
to as Papers I and II (Y. Suimazu [1955, 1956)), 
we studied chemical equilibrium conditions 
within the gravitational field and applied the 
results thus obtained to a detailed study of 
the chemical structure of the earth. However, 
Paper I was concerned only with the one- 

~phase and isothermal model of the earth, its 
constituents being FeO, MgO, SiO, and Fe. 
When Fe and O co-exist within a system, 
there occur two possibilities. Either Fe-O 
mixture phase is stable or FeO compound 
phase is stable. If it is the former, then 
under what physical (temperature, pressure) 
condition does it exist? If it is the latter, 
under what physical condition does this phase 
exist? In Paper II, we discussed this kind of 
phase stability and determined the pressure- 
temperature conditions for the phase stability. 
This condition is, however, a necessary condi- 
tion and not a sufficient condition, because the 
mol ratio of the components should vary with 
depth, as given in Paper I. Even if the tem- 
perature and pressure are sufficient to form 
FeO phase, gravitational separation effect due 


to density differences in costituents will drive 
the heavier Fe downward and the lighter O 
upward. Then the thicker FeO layer will not 
be formed because the mol ratio Fe:O is far 
from unity at both upper and deeper parts 
of the earth. Isolated O is predominant at 
the upper part and Fe at the lower. Follow- 
ing is a detailed study of the chemical phase 
equilibrium within self-gravitating planets 
based upon the results obtained in Papers I 
and II. 


2. Distribution of elements within self- 
gravitating planets 


Our simplified model of the planet has three 
independent components Fe-Si-O. Three 
chemical reactions written below are assumed 
to occur within the model planet. 


Fe+O2 FeO, (ik) 
Si+20 2 SiO. , (B) 
2FeO+ SiO, @ Fe.SiO, . (2) 


The compounds Fe,O; and Fe;0, are omitted 
from our consideration, because except at ythe 
surface layer of the earth they are apparently 
unstable (T. F. W. Barru [1948]). FeSiO; is 
also omitted since the thermochemical data 
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for FeSiO; reaction is not available and, fur- 
ther, the energy and density differences 
between FeSiO; and Fe.SiO, are neglected in 
our rough estimation. Reaction constants for 
(1)-(3) can be calculated as the functions of 
pressure and temperature using the results in 
Paper II. 

Chemical equilibrium conditions are written 
by the results in Paper I as follows: 


OG; OG; 
= — Sadr : dP 
iG te ee ae: ae 


+3(F") dny=0, (A) 
j On; PrD ts Mpeg 


where 

G:: molar Gisss’ free energy of 7z-th com- 
ponent, 

Si: molar entropy of z-th component, 

m;: molfraction of z-th component; 1=Fe, 
2=O0, 3=Si, 4=dependent ‘component 
(m=1—2,—N,—Ns), 

ry. distance from the center, 

P; pressure; 

T: temperature. 


Thermodynamical identities give 


(a = mig (5) 
Or Ba :L 5% 
and 
0Gi 
L =v, 6 

ee es . 

where 
m;: molar weight of z-th component, 


%: molar volume Of z-th component. 
Hydrostatic equilibrium condition gives 

OP 

Or 


Mean density p of the system is a function of 
molfraction of components and thus it depends 
upon 7 even if the incompressible planet is 
assumed. Nevertheless for simplicity, we 
put p=4.0. Transforming the independent 
variable 7 into 


=—p9. wa) 


a= 1—7 (a: radius of the planet), 


we get from (4)-(7) 


dinm _ __ Fs ga 
dz te eRT 


S, dT _AG,dn, 


RT dx RT dx ’ 
fins —(00.— ms) orn 
S, dT AG,dn, 1 4G, dns. 
TRI dz RTidg oR ee 
218 —(pvy— ms) Bo 
S; dT _14G, dn, 
RT dx 2RT dz’ 
where 4G, and 4G, are energy changes dur- 
ing the reactions from left to right in (1) and 
(2) respectively. Second terms of the right- 
hand sides of (8) vanish in the isothermal 
planet. We can neglect the second term for 
even 2,000°K (surface) 6,000°K (center), 
since it is only one-tenth of the first and 
third terms combined. It is noteworthy that 
dT/dx is the excess of the temperature 
gradient over the adiabatic gradient for the 
compressible planet. 
As shown in Paper II we have 


(8) 


AG=4AG"(T)+PadAv , (9) 
and 4v is given by 
4v=(va—vs) + , (10) 
OA 


under the reaction B+CT2 (See equation 
(14) of Paper ID. Values mi: and wv: of any 
material are given in Table I. Molar volume 
of O is calculated from ionic radius (T. F. W. 
Bartu [1948}). 


Table I. Molar weight m; and molar volume 2;. 
mi | Vi 
Fe 55.85 eo 
Si 28.10 12.06: 
O 16.00 6.93 
FeO 71.85 | 11.98 
| 
SiOz 60.19 | 22.64 
Fe,SiO, 203.8 | 47.61 


We get the unique solution of equations (8) 
if the molfraction 
mr ar:\ EY he dr: \"n dr (11) 


0 0 0 


for the whole planet is known. It is, how- 
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ever, obviously prohibitive to solve (8) 
numerically with boundary condition (11). 
We will begin integration by giving a suita- 
ble surface value. Trial and error method is 
used to satisfy 


s NM = 1 
t=1 
at any depth. 

If we do solve (8), we encounter the follow- 
ing serious difficulty. Assuming the surface 
temperature to be 2,000°K for the earth, 
equations (8) become at the surface 


dinn, dny 
——_-= + 102.6 .o0 —, 
cm +102.6+ 7 aa 
d\n No dn, dn, 
—=—44,4+7,80 + 15.35—— 
dx See dz * az 
d\n P13 dny 
Se 76.2-- 15.355 2 
“a (Oz Es (12) 


The first terms of the right-hand sides which 
give the effect of gravitational differentiation 
are considerably larger than the second terms 
which express the effect of chemical interac- 
tion. Then we see that Fe and O would 
separate completely if the structure of the 
present earth is formed under gs ~ 10° (=pre- 
sent value). The present value of gravity is 
too large! It appears that #, (surface) <1. 
The result of a rough calculation shows that 
m, (surface) must be less than exp (—16.5) 
stor 2 < 10¢°K. 


3. Variation of gravity during terrestrial 
evolution 

According to H. Brown and C. PaTTerson’s 
study (1948) on the composition of meteoritic 
proto-planet, the effect of gravitational dif- 
ferentiation is predominant in the deeper 
parts of the plant. The molfraction of O 
varies linearly from only 60 per cent at the 
top to 49 per cent at the bottom and variation 
of Fe is from 4 per cent to 15 per cent (see 


Fig. 2). From Fig. 2, we obtain 
d\n, d\n n, 
——— ~+1.1 a 6, 
dx 3 ; dx 
eH 0.16 (13) 
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at the surface. Putting (13) into (8) we get 
gsa/RT ~ 0.03 for the meteoritic proto-planet 
(gs: surface value of g). Then gsa must be 
0.02 X (Ys@)present earth for T=3,000°K, which 
is the temperature given in Brown’s paper. 
Brown deduced also that the chemical reac- 
tion in meteorites occurred at a high pressure 
of the order of 3x10° atmospheres. Since 


‘the present Mars has 3x10° atmospheres at 


the center, he came to the conclusion that 
the proto-planet would be of the same size 
as Mars. However, the present Mars has 
gsa~0.23 (9s@)present earth. From our point of 
view, this seems to be too large. 

In a rough calculation, gy can be replaced 
by the mean value gm. For a homogeneous 
planet gm is equal to 4 gs. Since gscca and 
the total mass M is Mca’, we can write 


gece (14) 


On the other hand. G. P. Kuiper (1949) 
discussed the variation of M/a during ter- 
restrial evolution. According to him, the 
lighter elements gradually escaped toward 
outer space after the terrestrial accretion has 
been completed, and M and a decreased 
gradually during terrestrial evolution. His 
study shows also that M/a was equal to 0.4x 
(M/d)present at the initial stage, decreased to 
the minimum value 1/40(V/@)present, and then 
increased to the present value. The escape 
of comparatively heavy elements such as Fe, 
Ni, Mg, Si stopped at the minimum M/a- 
stage and thereafter the variation of M/a 
was due mainly to the escape of chemically 
non-active gas such as Xe and Kr. The mini- 
mum value 1/40(M/a@)present agrees well with 
0.02 (9s@)present Obtained from (13). 

Are there other possible ways to avoid the 
serious difficulty already mentioned? Based 
upon the assumptions that chemical equili- 
brium is attained, we can assume 

(a) the temperature was much higher than 

LOSI) OF 
(b) hydrostatic equilibrium was not attain- 
ed at the stage of formation of chemi- 
cal structure. 
To satisfy assumption (a), the temperature 
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should be higher than 10°°K which seems to 
be too high. To satisfy assumption (b), we 
must assume that all chemical reactions oc- 
curred at the stage of large gsa and that 
afterwards the gravitational separation oc- 
curred at smaller gs@, which assumption seems 
to be rather artificial. 


ie] 


80 


60 


40 


20 


—s 
Fig. 1 (a) Distribution of elements within 2,000°K 
isothermal planet. 


Fig. 2. Distribution of elements 
within BROWN’s meteoritic proto- 
planet. (after H. BROWN and C. 
PATTERSON [1948]) 


4. Chemical constitution of the simplified 
earth’s model ; 

Solving equations (8) for the special case 
gsa/RT=0.03, we obtain Fig. 1 (a) and (b), 
which may be adopted tentatively. The dis- 
tribution of elements within Brown’s meteori- 
tic proto-planet is shown for comparison in 


——__-> 


Fig. 1 (b) Distribution of elements within 2,000°K 
(surface)— 6,000°K (center) planet. 
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Fig. 2. Substituting 2, 72, 23; thus obtained 
into (1)-(3), we get the complete structure of 
the model as shown in Fig. 3 (a) and (b). 


Fe2Si0a 


ee 
Chemical structure of 2,000°K 
isothermal planet. 


Fig. 3 (a) 


Chemical structure of 2,000°K 
(surface) > 6,000°K (center) planet. 


Fig. 3 (b) 


Fig. 3 (a) of the 2,000°K isothermal model 
shows that Fe and Si at any depth are of the 
oxide forms FeO and SiO,, and that the re- 
maining O is of isolated form. FeO and SiO, 
tend to combine to form Fe,SiO,. Since SiO, 
decreases with depth and FeO increases, Fe,- 
SiO, content shows maximum at about «= 
0.5. When isolated FeO and SiO, are at shal- 
lower or deeper depth than x«=0.5 respectively 
they will disappear. In Fig. 3 (b), i.e., the 
model with 2,000°K (surface) —6,000°K (cen- 
ter), FeO and SiO, decompose at x«=0.5 and 
x=0.8 respectively. Then a metal core will 
be formed of isolated Fe and Si. The mantle 
is composed of Fe,SiO,, remaining SiO, and 
isolated O. 

Since O content is too high in our simpli- 
fied model, it can be found at any depth. In 
the actual planet, however, isolated O will be 
of a combined form with elements other than 
Fe and Si, as for intance Al or Na. 


This work was carried out by the writer 
under receipt of Grant in Aid for Develop- 
mental Scientific Research of the Ministry of 
Education. 
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Temperature Distribution within the Mantle. 


By 
Michiyasu SHIMA 
Geophysical Institute, Faculty of Science, Kyoto University, Kyoto. 


Summary 


Temperature distributions within the mantle have been calculated by the modern 
theory of solid, on the assumptions that the thermal state is nearly stationary at present 
and there is only a conductive heat flow within the mantle. The main results obtained are 
as follows. Under the above assumptions, temperature distributions with Ye less than 2.25 
cannot be obtained in the B-layer (413-1000km). The temperature gradient there is 
nearly constant and increases gradually with depth in the D-layer. (1000-2898 km). 

Notations. g: density, oo: density at p=0 and T=0, p: pressure, 7: temperature, 
r: distance from the earth’s centre, g: acceleration of gravity, Kr: isothermal bulk 
modulus, Ks: adiabatic bulk modulus, «: coefficient of thermal expansion, Y¢: GRUNEISEN’s 
parameter, ¢: Ks/0, Vp: velocity of P-wave, vs: velocity of S-wave, A: mean atomic 
weight, N: LOSCHMIDT’s number, 4, u,m,m,l: LAME’s constants, K: thermal conductivity, 
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P: rate of generation of heat per unit volume, Cy: specific heat at constant volume, C>p: 


iN) 


specific heat at constant pressure, 0: 


There are two points of view concerning 
the thermal history of the earth. The first 
assumes that the earth started as a hot gas 
and has been cooling ever since. The second 
assumes that the earth started as a cold body 
and got warmer in the course of time. But 
if the earth had once entirely molten, its 
subsequent history would be the same as that 
of the earth with a hot origin. The point of 
view that the earth has been cooling through 
the molten state is considered correct by 
many geophysicists (B. GuTENBERG, 1951). The 
epoch immediately after the molten state was 
characterized by a fast cooling as the results 
of mixing process and radiation and the cool- 
ing velocity became very small after the 
mantle had solidified. Particularly, the cooling 
in the mantle and the deeper part is very 
slow (H. Jerrreys, 1952). Then the thermal 
state may be nearly stationary at present and 
thus it is assumed in the present calculations 
that there is only a conductive heat flow in 
the mantle and radioactive matters distributed 
uniformly there and the rates of heat gene- 
ration per unit volume are constant respec- 
tively in the homogeneous B- and D-layers. 

If the condition of hydrostatic equilibrium 


constant relating to crystal structure. 


is satisfied within the mantle, that is, dp= 
—godr, the variation in density in the homo- 
genous layer is expressed by the following 
equation (F. Bircu, 1952) 


dp _—/{ Op \- ap 00\ aT 
57 =( ap b dr Aa) dr 
et 9p aids 

a 


(1) 


By the thermodynamical relation between 
isothermal and adiabatic bulk moduli 


Fi il ol Flare 
Kon i 
the density variation (1) is written as follows 
ae a Stes eee = 
Hts toed gal de 
J0 
aa UL) (2) 
F; 0 


We shall divide the mantle into thin shells 
in the following numerical integrations, the 
outer and inner radii being ™ and 7%, and use 
the suffixes 1 and 2 to distinguish the physical 
quantities of the materials existing at 7% and 
7, (H. Mix, 1952). Then (2) may be approxi- 
mated by the following equation 
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_  20—p/o) 29 

7 (n—r)+ pal) | b+ $2" 
The ratio of densities o./o, can be derived 
from the definition of GRUNEISEN’s parameter, 


(3) 


( ape Ue +20s*)s (4) 
Or (Up 2+20s-%). ” 
where 
0 log vm 
re A log (Alo) ’ 
rae (N0.Y" (vp73-+-2u578)-1/3 = Coe ‘if . 
(5) 
Vp? =([A+2u+ f (164+ 20u—18/—42) 
+ f?(62.524+652—207/—55 —37)|/0 , 
(6) 


vs? =["+F (8A4+9+1.592+0.522) 
+ £7(152+27.54—271+1.5m+2.5n)]/o . 
(7) 
By (5), (6), (7), we have (Y. Suimazu, 1954) 


a os (30+Bif +B : 
where 
_(0/)?/8#— 1 
ue 9 ; 


B,=124+3(a/u—1)+ 1.572 +0.572 
723 


epee 1.592-+-0.5% 
a 


b 


’ 


B,=42.54+15(4/u—1) 2 7(1.522-++0.57) 
Ll 
(in the mantle) , 
(8) 


We use the power m-+-n=6 in the following 
expressions of the potential energy between 


=45.19 + 7(1.52+-0.592) 
Lu 


two atoms =“ ats B 


, yn 


(M. Suma, 1955). Then 02/0, are obtained 
from (4) and (8), taking rg as a parameter. 
The stationary and conductive heat flow is 
expressed as follows 
u e (nK 7 )=-P. 
TAT? dr 


We obtain by the integration of (9) 


, in deriving B, and B, 


(9) 


ai, =(% ) Ky _ (7? —70*) (10) 
b 


dr dr rk 3r2kK ’ 
where the suffix 4 indicates the physical quan- 
tities at lower boundaries of each layer. The 
heat conduction is considered non-metallic in 
the rocks which constitute the mantle. Thus, 
the lattice vibrational thermal conductivity A 
is expressed as follows (S. Mizusuima, 1954) 
@/3yV,C? 
ff ’ 
We put (11) in (10) and take (2) and (10) into 
consideration, then 
a. eae ae ie 
a PPP yo \ Ahy.ls 
EIS a) 
3r2dd23y pc? * 
On the other hand, from the following equa- 
tion 


K=0 d=p/a . (11) 


(12) 


La 1 aeee - _¢aA 
Cp C(1+Tare) ’ Cp 
we obtain 
cacel S (13) 
a T elo 
By putting (13) into (2) and (12), we have 
af) of Sten _ _ 7ndA 
mt oo 
“He or Tre) 
Cv Ta Q 


70°12 er? (Z)F =i A es os My 7° 
ry? \dr}]T, 370d7/3v,c? ~ 
(10) and (14) give on integrations 


T= eFn( | "0 eP ee vies ) 


’ 


r Talo 
Pir) = a (ate —are)dr. (16) 
is @ 
log T= log T ea Di Vylr*("a abr 
ai ar b AG r PER 


iz. R—7,3 


Se 39 \ Pd2/3y ,¢2 Te (17) 


where @ indicates the physical quantities at 
upper boundaries of each layer. @ and ¢ can 
be calculated from the observed values of vp» 
and ws and g is nearly constant within the 
mantle (10° cm sec-?). We eliminate P/éd and 
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A/c» from (15), (16) and (17) and calculate the stant in each homogeneous layer. The results 
temperature distribution, holding them con- are given in Table 1. 


Table 1. Absolute temperature within the mantle. 


st \e=2.25 ¥e= 225 ree | Yg=1.4 Ye=1.5 | 1a=2.0 
33 *1000 ; *1000 | 1000 *2000 *2500 *2000 *2500 | *2000  *2500 
100 JS 1142 1200 2102 2628 | 2190 2103) 9) 2369) 2898 
200 | 1270 1320 | 1400 2218 2762 2402 2935 2793 3340 
300 | 1380 1491 GOO ye 9 2843 2897 2646 3216 | 3230 3778 
413 1531 77S 1800 2470 3031 2920 3520 3700 4075 
| 2000 | 2596 S70) | 6201 3842 4259 4650 
2200 2730 3311 3537 4217 4802 5192 
2400 | 2887 3460 | 3903 4620 5396 5817 
2600 3068 3616 | 4314 5098 6000 6408 
2800 | 3280 3810 4720 5545 6595 7028 


* Assumed values. 


The values of JT are plotted in Fig. 1 against ascertained in it. The temperature distribu- 
the depth. tions in the B-layer with y« less than 2.25 are 
not given in Table 1, because we cannot 
obtain them keeping P/d and A/c, constant in 
it. The temperature gradients obtained are 
7000 a larger than the adiabatic one, and then the 

; heat may be carried away by conduction. 
They are nearly constant in the B-layer and 
gradually increase with depth in the D-layer. 


1G 
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Thermal Process Near the Earth’s Crust. 


By 


Keiiti Axt 


Geophysical Institute, Faculty of Science, Tokyo University 


Summary 


In connection with the earthquake energy problem, the thermal processes which can 
result in the steady accumulation of necessary energy somewhere close to the bottom of 
the earth’s crust has been investigated. With the aid of experimental and theoretical 
laws concerning thermal convections, some informations about the physical state at the 


place have been obtained. 


Introduction 


It is widely believed that an earthquake is 
a phenomenon of release of strain energy 
which has been stored up in the earth’s crust. 
The estimated energy of the largest possible 
earthquake exceeds 1074 erg. It is a very 
interesting and important problem how and 
why such a tremendous amount of energy 
can be accumulated and stored up within the 
earth’s crust. In this paper the writer is 
going to investigate thermal processes which 
are possibly taking place near the earth’s 
crust and which he believes are closely con- 
nected to the present problem. 

As to the fundamental cause of an earth- 
quake, there are at least two important hy- 
potheses. In one of the two, the stress genera- 
ted within the earth’s crust by subcrustal con- 
vection currents plays an important role and 
this view is supported by D. Griaes and Vening 
Mernesz, while in the other the change in 
internal pressure accompanying phase changes 
or general chemical reactions of involved 
material is essential and this latter view is 
supported by J. Jory and R. A. Dany. Recently 
T. Matuzawa (1953) has investigated the matter 
from the latter point of view and made 
quantitative discussions how a _ solid-liquid 
phase transformation can supply strain energy 
to the crust and make an earthquake occur. 
According to him if a subcrustal material 
having a volume of 500 km? melts being ac- 
companied by about 10 % increase of its own 
volume a pressure increment of 300 bars at 


the place will be generated. Thus the strain 
energy of 8x10%4 erg will’ be stored in the 
crust. The work done by the elevation of 
the crust through the gravity field is 4x10” 
erg. On the other hand the thermal energy 
consumed as the latent heat is estimated to 
be 9x10?" erg. In his paper Prof. Maruzawa 
was chiefly concerned with the mechanical 
aspect of the problem, and the question how 
the thermal energy can be supplied to the 
place concerned still remained open. Before 
studying quantitatively about the possibility of 
the supply of heat energy to the field of 
earthquake, it is important to select such 
processes that will satisfy a certain condition 
among many processes which may be assumed 
to take place under the crust. As to this 
condition, that for stationary increase of heat 
storage will be reasonably taken. 

The conversion of thermal energy into 
mechanical energy within a system of material 
usually takes place through the change in 
volume of the system; for instance, even the 
thermal convection current in the gravity 
field is maintained by the work done by 
volume change. The volume change of a 
part of material is a result of the accumu- 
lation of heat into that part. If there is no 
such accumulation, the conversion into me- 
chanical energy can never take place however 
abundant the heat flow through the material 
may be. In the following discussions we shall 
confine ourselves to investigate those processes 
in which heat flows steadily upward from a 
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high temperature source assumed to exist at 
greater depths of the earth and the heat is 
accumulated somewhere close to the crust 
steadily. 


Condition for the steady accumulation of 
heat 


The condition for the steady accumulation 
of heat into a part of material is satisfied 
only in the following two cases. Putting the 
velocity of flow at the part V, 

> ¥2<0 

2. or in the case of V=0, denoting the 
number of components in the system by 8, 
and the number of phases by a, 


a=B+1. 
We can verify the first of these statements 


easily. The accumulation of heat into a part 
of material from the outside per unit time is 


Rk | p>Tdav, (Ll) 
where k& is the thermal conductivity, The 
equation of heat transfer is 

OT .(V -grad)T=ap'T, a= : >» ad 
ot oc 


where c is the specific heat and o the density. 
Therefore even in the steady state OF eh. if 


V0, the heat accumulation can take place. 
This case corresponds to the process required 
for convection theory of the earthquake cause. 

In the next place, as to the second of the 
statements the accumulation will appear to be 
es =. But 


transfer is correctly 


impossible in the case V=0, and 
the equation of heat 
written as 
GS oF 
T—= 
ot 
where S is the entropy per unit mass. On 
the other hand there is the following relation, 


se Ox  (OS\ meer 
6: OT) 2g Ge ee dt 


+32( eS ) 0M, 
~"~X\OM,/e.7. Ot 


where FP is the pressure and M;, is the number 


yr ¢3) 


(4) 


of mol of the 2-th component of the z-th phase 
of the system. Therefore it may happen that 
even if temperature and pressure are kept 
constant, the phase transformation enables the 
heat accumulation into the system. The com- 
position of a phase can be determined by its 


concentration c,, and we have 


c= Mz/> M,=M,/M! (5) 
OMi dc, . « OM! 
fiona (6) 


Assuming the system as a closed one, we get 


0 F 
ap (as Mn)=0 (n=1, 2, ---- B) (7) 


We shall consider the case of a<8+2, for 
the case of a> §-+2 cannot exist by the phase 
rule. At first, in the case of a=8+2, the 
pressure, temperature and concentration should 
take the thermodynamically determined values 
for the system. But it seems unnecessary to 
consider such a strict condition at the interior 
of the earth. 

In the case of @a=8+1, if the value of 
pressure is once determined, the temperature 
and concentration are also determined at the 
same time. Therefore if P is kept constant 
by some external condition, the temperature 
and concentration of the system are automati- 
cally kept constant as long as the condition 
a=8-+1 is fulfilled, and we have from the 
equation (4) and (6) the following equation 


0S se a) oT +( OS OP 
+>'> ( = ) pen 
=e GME fre Ot 
OS 


OM 
yl é ‘ 
ss 5 Pa omMi jee ot 


SPP ca “alae 
& le OMi,/e.r. ") Ot 


On the other hand, we have from the equation 
(8) 


(9) 


. om 
G . 


=Q “red 98 Sry 
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These equations (9) and (10) together represent 

a set of linear equations about the variables 
OM 
ot 


The number of equations is 8+1, and it is 
equal to the number of variables in the case 
of w=8+1, therefore we can have not vanish- 


aM as 
eT Then. the 


G=1,2,+---, a), 


at and also ate oT 


steady accumulation of heat is possible to 
occur in the case of a=8+1. 

Finally in the case of a< 8, if the stationary 
condition is established i.e. the pressure, 
temperature and concentration are kept con- 
stant, a set of linear equations can be derived 
as above. In this case, however, the number 
of equations 8+1 exceeds the number of 
variables a. And it can be readily seen that 
the simultaneous equations are satisfied only 


Thus in 


ing - 


; OS : 
when all the M,’s and vanish. 


Ot 

this case there should be no steady accumu- 
lation of heat into the system, and the veri- 
fication of the whole statements is completed. 

For example, in the case of one single 
component, the steady accumulation of heat 
can take place only in the system in which 
two phases exist at the same time. If the 
temperature of material at the bottom of the 
earth’s crust is at the melting point, then it 
should have a constant value determined by 
the value of the pressure at the place, and 
consequently the quantity of heat flow through 
the crust should be determined. Therefore if 
the heat flow from below into the place exceeds 
the above determined value, the extra heat 
will be accumulated at the place and the con- 
siderable part of extra heat will be converted 
into mechanical energy. Hereafter we shail 
call the substance at the place “working 
substance.” 
Heat supply for the phase change 

The heat flow to the earth’s surface from 
below is estimated as 1.3 10-° cal. cm™ sec7! 
by the observation at the surface. On the 
other hand, if an earthquake of energy of 
0(1025) erg occurs in a circular area of radius 
50km once in a hundred years, Matuzawa’s 


calculations have shown that heat energy of 
O(10*") erg should be accumulated into the 
working substance which is assumed at the 
bottom of the crust of the area during the 
period ; this is equivalent to the heat flow of 
OQ(10-*) cal. cm-? sec-!. This amounts to a 
hundred times of the quantity of heat flow 


tH. 


Heat flow through 
the working substance. 


Bigg 1h 


through the crust. We shall consider in what 
way such a quantity of heat can be supplied. 

A) radioactive energy; At first we shall 
investigate the possibility that the necessary 
heat can be generated by radioactive substances 
contained in the working substance itself. 
Table 1 shows the radioactive heat generation 
in some igneous rocks (Bircu, F: 1942). The 


Table 1. Radioactive heat generation of igneous 
rocks. (after R.D. EVANS and C. GOODMAN) 


Rock Type Heat Generation 


Acidic Igneous Rock | 1-7 10=8 caller. ae sece 
Basic Igneous Rock | 0.6% 10=8 
Ultrabasic Igneous Rock 0.3 10-8 


volume of the working substance is assumed 
as 500km? in Maruzawa’s calculation, so 
that in a hundred years this may generate 
the thermal energy of O(10") erg and this is 
very little compared to the necessary energy 
of O(102") erg. Therefore if there is no extra- 
ordinary concentration of radioactive sub- 
stance, it is impossible to maintain the process 
by the heat generated inside the working 
substance. Thus the necessary heat should 
be transfered from the outside by means of 
thermal conduction or convection. 

B) thermal conduction; At the earth’s 
surface the heat flow from below is observed 
to have the mean value of 1.3x10-° cal. cm~ 
sec-!. If we take as the value of thermal 
conductivity of the crustal substance 0.006 
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cal sec“! cm-! deg-!, a reasonable value, the 
temperature gradient throughout the crust 
may be assumed to be 20°C km-!. There 
are two possibilities to enable the accumu- 
lation of necessary heat only through the 
thermal conduction into the working sub- 
stance. If the thermal conductivity beneath 
the working substance has the same value as 
that in the crust, the temperature gradient 
there should be a hundred times of that in 
the crust. On the other hand, if the temper- 
ature gradient beneath the working substance 
is the same as that throughout the crust, the 
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Fig. 2. Thermal conductivity of Basalt 
(after POOLE). 


thermal conductivity there should be a hundred 
times of that in the crust. 

The first possibility is difficult to conceive, 
for the temperature of a hypothetical heat 
source should exceeds 10,000°C, if it is taken 
to be located at the place 5km apart from 
the working substance. 

Then the second possibility that the thermal 
conductivity at the place beneath the bottom 
of the crust can be a hundred times of that 
throughout the crust should be examined. 
The temperature dependence of the thermal 
conductivity of basalt and granite has been 
studied by Pootr (1921), the results being 
shown in Fig. 2 and 3. The increase in 
thermal conductivity with the temperature 
increase is observed only in the case of basalt 


at lower temperatures. Generally the thermal 
conductivity of material decreases with the 
temperature increase except in the case of 
graphite. hematite, quartz glass and asbest. 
The effect of pressure is to increase the 
thermal conductivity. The pressure dependene 
of the thermal conductivity of basalt has been 
given by P. W. Bripeman (1949) as 


K=0.00404+0.000019 (P/1000) 
(P in bars. K in cal. sec-!. cm~! deg-?) 
=0.00414+0.0000089 (P/1000) 
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Fig. 3. Thermal conductivity of Granite 
(after POOLE). 


The increment is only 2.1~4.7% by the 
pressure of 10,000 bars which corresponds to 
the depth about 35 km below the surface, and 
at most 50% by the pressure at the depth 
of 300km if the equation holds to such a 
high pressure. Thus it is not probable that 
the thermal conductivity of substance under 
the crust will be notably high owing the effect 
of high temperature and pressure at the place. 
The change of material itself with the depth 
may increase the conductivity at the place, 
but as shown in Fig. 4 drawn based on the 
data given in F. Bircn’s table, the difference 
of the values among these rocks at lower 
temperatures gradually diminishes with the 
increase of temperature. It seems reasonable 
therefore to conclude that the value of thermal 
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conductivity under the crust is of the same 
order as that throughout ‘the crust, and the 
possibility of the required heat accumulation 
into the working substance only through the 
conduction should be rejected. 

C) thermal convection: The problem of 
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Fig. 4. Thermal conductivity of some 
rocks and minerals (from BIRCH’s table) 


heat transfer by means of the thermal con- 
vection has the long history of theoretical 
and experimental studies (FiscueNpDER, M., and 

Saunpers, O. A.: 1950; Ecxart, E.: 1949). The 
~ experimental study stands on the basis of 
dimensional analysis made by W. Nusseu7, and 
aims to establish the relation among the follow- 


Table 2. 
shape | l 
vertical cylinder height 
vertical plate | height 
vertical plate or vertical cylin- : 
der of large diameter | height 
horizontal plate breadth 
horizontal cylinder diameter 


horizontal or vertical cylinder | diameter 


two horizontal plates distance/2 | 


two vertical plates 


ing nondimensional numbers: 


Al 

ko 

Here / is the length representing the dimension 
of convection, @ is the temperature difference 
between the fluid and the boundary, k is the 
thermal conductivity and H is the heat flow 
carried by the convection current. Nz re- 
presents the ratio of convectional heat flow to 
conductional heat flow. 


i) NussELT number: Na= 


ii) ReryYNoLDS number: Re—Uet 


LL 


ili) PRANDTL number: iets 


in which ¢ is the specific heat and y is the 
viscosity. This is the ratio of the viscous 
diffusion coefficient to the thermal diffusion 
coefficient. 


Gye lee 
we 

in which a is the thermal expansion coefficient 
and yg is the gravitational acceleration. This 
number has the same character as that of the 
REYNOLDS number in the case of fluid motion 
forced by the buoyancy. We assume that the 
subcrustal convection current is not of the 
nature of forced convection which is supported 
by some external stress field but of one of 
free convection supported by the buoyancy 
generated by the thermal expansion of sub- 
stance. In this case we need not take into ac- 
count the ReyNoLtps number, and a relation can 
be deduced among Nz, Pr and Gy. It has been 
derived by theoretical considerations, though 


iv) GrasHor number: 


authorities 


c | c! | 
Onno | 0.067 TOULKIAN, HAWKINS, JACOB 
0.64 _ Boscu, McADAMS 
0.56 0.12 | SAUNDER 
0.54 0.14 | SAUNDER 
0.53 BoscH, McADAMS 
0.47 0.10 SAUNDER 
0.44 0.19 MULL, REIHER 


| L\"/o L ye aa 
jedi SM) So me | MULL, REIHER 
jdistance/2 | 0.41 ey | 0.18 ie | 


1) EEE Eee eee 
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there are some objections, that Nw is the 
function of only PrxGr. We have the follow- 
ing formula in the case of laminar flow, 


Nu=c(Pr-Gr)'/* (11) 
and in the case of turbulent flow 
Nu=c (Pr-Gr)? (12) 


The experimental studies have determined c 
and c’ in many cases of boundary of different 
geometrical shapes. Table 2 is a summary of 
some of these experiments, in which the 
material used are air, water, oil and ethylen 
glicol. 

Although we do not know about the geometry 


0 
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Fig. 5. Relation between the NUSSELT number 
and the RAYLEIGH number for several geo- 
metrical shapes. 


of the heat transfer actually taking place 
under the crust, it is reasonable to assume 
that the actual case is included in the range 
of variation of these examples. If the neces- 
sary heat accumulation into the working sub- 
stance is done only through the convection 
current, and if the thermal conductivity and 
the gradient under the crust are assumed to 
be of the same order as that throughout the 
crust, the Nussett number Wz, i.e. the ratio 
of convectional heat flow to conductional heat 
flow, should be taken 100. The value of Prx Gr 
corresponding to Nussetr number 100 is in 
the range of 108°~10°° as shown in Fig. 5 


which is drawn using the data of Table 2 in 
the case of laminar flow. 

Now we shall investigate how much can be 
said as to the physical state under the crust 
from the requirement that Prx Gr should take 
the value of 10°-°~10°:° there. 
cag o*l* 

ku 


PrxGr= (13) 
has the name of “RayLeIGH number”, and 
determines the condition of generation of 
Benarp’s cell. Among the physical constants 
at the right hand of the above equation, the 
specific heat c and the thermal expansion 
coefficient a have somewhat similar nature; 
both have zero value at the absolute zero 
temperature, increase with the temperature, 
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Fig. 6. Temperature dependence of 
heat capacity of some rocks (from 
BIRCH’s table). 


and take extraordinary large value at and near 
the melting point. On the other hand, the 
viscosity # and the thermal conductivity k 
are both related to the transport phenomena; 
the former is concerned with the transport of 
momentum and the latter with that of heat. 
As the similar quantities are connected by 
the multiplication in the definition, the Ray- 
LEIGH number is expected to be very sensible 
to the change of the state. 

Among the physical constants involved in 
the RayLEIgH number, we can estimate from 
laboratory experiments the values of specific 
heat c, thermal expansion coefficient a, density 
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o, and thermal conductivity k of the substance 
under the crust. , 

Therefore, if we put the temperature 
gradient 6/7 to be the same as that in the 
crust, we can obtain a relation between the 
two unpredictable parameters, i.e. the viscosity 
# and the linear dimension of convection 
current 7, from the requirement that the 
RayLEIGH number should be within the range 
Of 108-210? 

The specific heat has nearly constant value 
above the Dresyr temperature except in the 
neighbourhood of melting point. Fig. 6 shows 


10% pecree! 
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Fig. 7. Temperature dependence of 
thermal expansion of Diopside. 
(from BIRCH’s table) 


Table 3. Thermal expansion of rocks in 


solid and liquid state. 


« (at ordinary 


Rock Type temp.) a (liquid) 
Basalt | 1.5x10-5deg-1 | 8.2x10-5deg-! 
Diorite | 2.1 14.0 


the temperature dependence of the specific 
heat of basalt, diabase and granite. And the 
average value of the specific heat of rock at 
high temperatures, excluding the neighbour- 
hood of the melting point, is estimated as 1.3 
joules gr-?. deg.~1. 

As to the coefficient of thermal expansion, 
the temperature dependence in the case of 
diopside is shown in Fig. 7. The value at 
and near the melting point is extrordinary 
large as in the case of the specific heat, and 
the value in the liquid state is somewhat larger 
than that in the solid state. Table 3 shows 
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the average value of thermal expansion of 
some rocks at ordinary temperatures and in 
the liquid state, though the samples used are 
different between the solid state average and 
the liquid state average. 

The thermal expansion coefficient shows a 
considerable change by the effect of pressure. 
The figures in Table 4 is derived from the 
result of experiments made by P. W. Bripaman 
in the case of alkali metals, 4v is the change 
of volume and da is the change of thermal 
expansion coefficient by the effect of pressure. 
The ratio 4a/a) takes almost same value for 


Table 4. Change of thermal expansion by 
the change of volume for Alkali metals 
(estimated from BRIDGMAN’s data) 


— Av/V Aaa (Li) | Aa/ao (Na) | Aalag (EK) 
0.05 | 0.24 0.20 0.17 
0.10 0.42 0.37 0.35 


a fixed value of 4v/v. among the three alkali 
metals. If it is allowed to extrapolate these 
figures to the substance in question in the 
corresponding state, the value at the earth’s 
surface should be reduced by half by the 
pressure at the depth of 300km. Under these 
circumstances it seems reasonable to take as 
the value of thermal expansion of the substance 
under the crust 5x 10-° deg.-! unless the state 
is at the melting point. 

The value of the thermal conductivity, of 
which the temperature and pressure de- 
pendences are refered to before, is assumed 
to be 2x10-? joules sec.-! cm-} deg.-1, and 
also the value of the thermal gradient 20°C 
km-!, and the density 2.6. 

Using the above numerical values, we can 
obtain the following relation by the require- 
ment that the RayteicH number should be 
LOR P10? 2; 


Pee BAK OB LO (14) 


Here the units of 7 and of w are km and 
poise respectively. To enable the accumu- 
lation of the necessary quantity of heat into 
the working substance, the above relation 
should be satisfied between the dimension of 
the convection current and the viscosity of 
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material under the crust. The diagram of 
/4--1=10-* is drawn in Fig. 8. 

Thus the necessary heat energy for a large 
earthquake can be accumulated to the working 
substance, if the relation expressed in the 
formula (14) holds at the place under the 
crust between the linear dimension / of con- 
vection current and viscosity 4, for instance, 
7=10km and ~=10" poise, or 7=100km and 
m= 10" poise: 


log, + (poise) 


10 $0 100 


Fig. 8. Relation between the linear 
dimension of convection current and 
the viscosity of material for the 
NUSSELT number 100. 
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Fig. 9. Assumed velocity profile. 


We have not been able to say anything 
about the velocity of convection current in 
the foregoing discussions, but if we assume 
a simple model as to the pattern of current, 
we can estimate the value of the velocity 
and furthermore of stress generated at the 
bottom of the crust by the current from the 
requirement that the NusseLtt number should 
be 100. 


Estimation of velocity of convection current 


Let us consider a simple model of heat 
transfer from fluid to solid wall through 


stationary convection current as shown in Fig. 
9, in which z-axis is the boundary between 
the wall and the fluid, and the fluid at high 
temperature ascends along y-axis from below, 
and flows turning to the both sides along z- 
axis. The equation of heat transfer is 


DT k 
=ap? =— 15) 
De Ds, 4a au (15) 
if we assume near the wall 
oT oT oT eT 
” Ox * ” ay 2G Oy? > Ox, 
then the equation (15) reduces to 
aT 2 ork 
u ae a aye” (16) 


Developing z in Tayior’s series near the wall, 
and taking only the first term, we have 


Ou : 
u=( ay ta BY f 


taking this into account we obtain the follow- 
ing formula from the equation (16) 

OT 

Ser ; iy 
- Oy? ( ) 
Changing the variables 2 and y into w ac- 
cording to the relation w=ya-/3, we have 
the following equation with a single inde- 
pendent variable, 


Sab ras lars see 
dw? * 3a “ dw Ns (18) 
The solution is given by 
< = B 8 
T=K,\ e% dotk,. (19) 
0 
If the temperature of the wall is %, and 


that of the ascending fluid 4, the boundary 
conditions are given by 


y=0, o=0, T=h, 
and 
oar = cok seit 
Then we have 
K,=th, 
and 


Katt / | e Si dw 
0 
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= (4 I (t:t) ( eB he 


ad 


The temperature gradient at the wall is 
1\1/3 
es (3) 
Oy al 1 
i Rs i 
© 


The heat flow in the direction perpendicular 
to the wall per length Z of the wall is 


o=-fa(2) 


1 Nive 
3 
z 45) 


B\iW¥3 
.-4)(7) ale (20) 


B 
a 


7, ae 
mt—t)( ) Ts, (21) 


and the NusseELT number can be expressed as 
Q-1 7 (2 1/3 : 

=— =0.8x — 
— LR(t,—tp) OL =) pe) 


by which we can estimate the value of B, 
i.e. the velocity gradient at the wall, for a 
given NussELT number, and 
Bey. L Ga) 
Zé \ 0.8 
The value of a is obtained by the foregoing 


Nu 


(23) 


estimation as 0.6x10~ cm?.sec.-! Let Z be 
equal to J, and Nz 100, we have 
B=12X10-* x1 sec? (24) 


in which the unit of Z is km. 

The tangential stress at the wall generated 
by the current is 

t=p~B=O0(105)x/? dyne cm-?. (25) 

For the estimation of the absolute value of 
current velocity by the above simple method 
we need to know about the type of distri- 
bution of velocity, and we assume 


y y\? 
~ 122. i 
u ut ( a) 


This represents the distribution of flow in 
which the velocity in the direction parallel to 
the z-axis is zero at the wall and y=/, and 
some authors have used this type of velocity 


(26) 


distribution to fit to their results of experiment 
about free convection current. In this case 
4 

Bie 
27 
so that we can estimate the expected maximum 
velocity under the crust as 


the maximum velocity Umax is equal to 


Umax =5 X 108 x 7-1 (m year-) (27) 
10*dyne/cm 
10 
A 
6 
4 
2 
km 
10 50 189 
Fig. 10. Relation between the linear 


dimension of convection current and 
the stress generated in the case of 
NUSSELT number 100. 


m/year 
1000 
500 
its ,9) 
ee 
10 50 100 km 
Fig. 11. Relation between the linear 


dimension of convection current and 
its maximum velocity in the case 
of NUSSELT number 100. 


in which the units of Umax and / are m year~! 
and km respectively. The relations (27) and 
(25) tell us that if the thermal convection 
current transfers and supplies the required 
energy to an earthquake, the velocity of the 
current and the tangential stress generated 


62 Keiiti AKI 


by it at the bottom of the crust are de- 
termined as the function of the linear di- 
mension of the current, as shown in Figs. 10 
ghavel: JUk 


Résumé 


The possible thermal processes near the 
earth’s crust are studied with special regard 
to a process presented by Prof. Matuzawa in 
his recent paper titled “Feldtheorie der Erd- 
beben.”, in which a solid-liquid transformation 
plays an important part. Before the discussion 
of this particular process, the writer gives 
the necessary condition for the stationary 
conversion of thermal energy into mechanical 
energy inside a system of material through 
which a stationary heat flow exists. It is 
shown that this condition is satisfied near the 
crust in two and only two cases, of which in 
the first a mass flow exists in the system, 
and in the second the relation w=8+1 exists 
between the number a@ of phases and the 
number 8 of components of the system. The 
first case corresponds to the convection theo- 
ries of the cause of earthquake and the second 
corresponds to the phase transformation theo- 
ries, to which Prof. Matuzawa’s theory be- 
longs. The necessary heat supply for a great 
earthquake was shown quantitatively in his 
paper. It is empbasized that such an amount 
of heat energy can only be supplied by the 
convection current under the crust. And the 
writer obtains by the aid of experimental and 
theoretical laws of heat transfer by a free 
convection the relation between the scale of 
convection and the viscosity of material under 
the crust required for the effective accumu- 
lation of necessary energy. If such a con- 
vection has a linear dimension of 100 km, the 


‘viscosity of material must be 10" poise, its 


maximum current velocity is estimated to be 
50 m/year and it exerts a tangential stress of 
the order of 10% dyne/cm? along the lower 
surface of the crust. 

Those results that the subcrustal current 
which is believed to account for the mountain 
formation may accompany energy accumu- 
lation by phase transformation, that the phase 
transformation which explains well the cause 
of a great earthquake needs the heat energy 
supplied by the subcrustal current and that the 
possible stationary conversion of the thermal 
energy into mechanical energy must be con- 
fined to the above two cases, may lead us to 
the conclusion that those two processes com- 
bined play a very important part in the 
physical phenomena near the earth’s crust. 
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Summary 


Various facts appear to suggest that one continuous field of mechanical stress deve- 
loped in the earth’s crust has a certain upper limit for its voluminal extent. The ultimate 
mechanical stress energy that can be stored up in this whole volume until a break-down 
takes place in it may be identified with the energy of the largest possible earthquake. 
The energy deduced on this hypothesis agrees well with those of the actual largest earth- 
quakes. The area A in which aftershocks occur in association with a major earthquake 
has been found by UTSU and SEKI regularly to increase with the magnitude M of that 
main shock. This relation, when combined with the magnitude-energy relation due to 
GUTENBERG and RICHTER, yields a formula 


E=6 x 102 x Al's, 


The numerical values of the coefficient and of the exponent of A in this formula can be 
well explained by the hypothesis stated above regarding the spatial distribution of the 
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stress energy within the earth’s crust. 


The amount of energy which is sent out 
from the ‘‘origin’’ of a large earthquake in 
the form of elastic waves is so enormous that 
it will be difficult to conceive of this much 
amount having been stored up within a small 
confined volume of the earth’s crust until the 
outbreak takes place”. Since the material of 
“which the earth’s crust is made up has the 
limit of strength which is finite, a huge 
volume of the earth’s crust must be needed 
for this much energy to be stored up in it in 
order that at no part within the volume the 
stress should exceed this limit of strength. 
Let us call this volume the ‘‘earthquake 
volume’’. As will be shown later, the stress 
energy that can be stored up within a unit 
volume of rock is about 5x10 ergs. 

On the other hand, according to our ex- 
perience throughout the long history of earth- 


1) Seismological Laboratory, Division of Geolo- 
gical Sciences, California Institute of Technology, 
Pasadena, Contribution No. 794. 

2) Discussions will be confined to shallow shocks 


only 


quake studies, there is likely a certain upper 
limit of the energy of one single shock, 
however large it may be. The largest 
magnitude of an earthquake listed in GuTEN- 
BERG and Ricuter’s book ‘‘ Seismicity of the 
Earth’’ is 8.6. Even in Japan which is noted 
for its high seismic activity and also for its 
richness in dependable historical records, no 
description can be found of an earthquake 
which was so large as to work damages over 
more than one third or one fourth of the 
whole area of the country. This may be 
taken to be another indication that the earth- 
quake energy or the earthquake volume has 
also a certain upper limit. 

In an article published in 1940, the present 
writer (Tsusot, 1940) proposed to calculate the 
energy & of a largest possible earthquake. 
First of all, he has assumed that the energy 
FE is given by 


: ex? V, Ge) 


2 
where the notations used have the following 
meanings: 


k= 
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e: effective elastic constant of the 
crustal material, 

x: ultimate strain of the crustal 
material, 

V: earthquake volume. 


The expression (1) amounts to mean that the 
material within the earthquake volume will 
be uniformly strained until this uniform strain 
reaches the value x everywhere at the same 
time within the volume, when, at the time of 
earthquake occurrence, the stored energy E 
is sent out from the whole volume by some 
kind of mechanism or other. Evidently this 
uniform strain hypothesis cannot be strictly 
true, but the essential point of it which the 
writer would like to emphasize is that 2 is 
assumed not to differ notably from one part 
to another within the volume. To take z 
everywhere uniform will be understood to be no 
more than the first approximation in the for- 
mulation of the hypothesis. According to 
the writer’s view presented in the article 
mentioned, vertically the volume cannot 
exceed the thickness d of the earth’s crust, 
nor horizontally three times the thickness, 
3d. This value of the factor 3 has been 
derived from the study of local versus re- 
gional isostasy. (Tsusor1, 1939). It is related 
to the ratio 
extent of regionality of isostasy 
thickness of the earth’s crust ’ 


where the extent of regionality means, after 
Vening Metnesz, the horizontal extent of the 
largest topography that can be supported by 
the strength of the earth’s crust, that is by 
a mechanism other than isostasy,. 

There is also another information in favour 
of the finite horizontal extent of the earth- 
quake volume of a largest possible shock. 
In 1949, the writer (Tsusot, 1949) investigated 
how the correlation coefficient between the 
year-to-year numbers of earthquakes (M>5.5) 
in various compartments of suitable size 
taken along Japanese islands decreases with 
the distance a between the two compartments 
compared. The compartments are bounded 
by straight lines which are drawn across the 
islands perpendicularly to their general trend. 


The correlation coefficients was found to 
decrease very regularly with the distance 
until it becomes nearly zero at a@=150 km. 
The occurrences of these earthquakes. may 
be understood to be the indications of the ac- 
cumulation of energy ia a volume V of the 
earth’s crust, of which the earthquake volume 
v for each one of the shocks is a small part. 
That the correlation coefficient becomes zero 
at @=150km would mean that the volume V 
cannot extent beyond that length horizontally. 

From the above considerations, & may be 
written as 


E=5 ex?d-3d-3d= 3 éx7d*. (2) 


The value of e may be taken to be 
5x10%~10" C.G.S.2X The value of zx has 
been known to be 10-4~2x10-* from the 
studies of deformations of the earth’s crust 
caused by large earthquakes, particularly by 
those which were associated with seismic 
faults (Tsupor, 1933). Thus 1/2 ez? which 
is the maximum energy to be stored up with- 
in a unit volume of rock will be 1/2 (5x10!) 
x 10-8 ~1/2 (10#2)« (4x 10-*) or 3%10°~2x10t, 
The thickness d of the earth’s crust may be 
taken to be 4x10°~5x10°. Putting these 
numerical values into the expression (1) for 2, 
the values 1.4x107*~2.3x1025 are obtained. 
If we take e=5 X10", z=10-*; d=5.x10*, the 
value of # will be 3x10”. 

From the nature of the problem, no great 
accuracy can be claimed of these values. 
They may well be uncertain by a factor of 
even as much as 5. 

Now on the other hand, according to the 
1942 formula of B. GuTENBERG and C. RIcHER, 
(1942), the energy of an earthquake is related 
to its magnitude M as follows: 


log E=1.8 M+12. (3) 


If we put M=8.6 into this formula, which 
is the largest magnitude listed in GuTENnBERG 
and RicuTEr’s book ‘‘Seismicity of the Earth,”’ 


E=3 x10" 


is obtained which value is far larger than 


9% 


3) Hereafter, all the numerical quantities will 
be expressed in C.G.S. system. 
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any of the values 1.4x 14%~2.3 x 10% estimat- 
ed by the present writer. Thus GuTENBERG 
and Ricuter have naturally stated in their 
book as follows: 
teteee This (Tsusor’s estimation)” is ap- 
preciably less than the largest energies cal- 
culated from the revised equation (2) (log E 
=1.8M+12) for observed earthquakes since 
1904, which appear to reach 102” ergs------ ae 

According to the Progress Report, Seismo- 
logical Laboratory, California Institute of 
Technology, 1955 (Bentorr, GUTENBERG, PRESS 
and RicuTer, 1956), however, the two inves- 
tigators have recently revised the formula 
into 

log £=1.5 M+11.8. 

If we put M=8.6 into this new formula, 


E=5 x 1024 


(4) 


is obtained. This value agrees well with the 
writer’s value 3x10". The agreement may 
be said to be very good considering the nature 
of the problem and the uncertainties of the 
numerical values involved in the calculations. 

Although the present writer feels himself 
much satisfield by this good agreement, he 
is furthermore going to see if there is any 
other additional information which will sup- 
port his hypothesis concerning the spatial 
extent of the earthquake energy field. 

In 1954, T. Ursu and A. Sexi published an 
interesting article in which they studied the 
relation between the horizontal area A in 
which aftershocks of a large earthquake 
take place and the magnitude M of the main 
shock. Using as data 30 large earthquakes 
-~ which had taken place in Japan, they obtained 
a formula as follows. 


log A=M-+6. 
Let us call A the ‘‘aftershock area.’’ 


If, for instance, we put M=8 
formula, 


(5) 
in this 


A=10"= (100 km)? 


is obtained. Combining this formula of Ustu 
and Sexi with that of GuTeNBERG and RICHTER 
(4), M can be eliminated with the result 


4) Parentheses by the present writer. 


log E—1.5 log A=2.8 
or E=6x10?x Al, (6) 


It has already been stated that according to 
the present writer’s view, E is given by 


E= (2) 


If the aftershock area A is the earthquake 
volume projected on the earth’s surface, 


which assumption appears to be a reasonable 
one, A is equal to 9d?, 


9 
2 3 
9 et a, 


it 
Sie 1*5 
or a= 27 Als, 


Thus we can write 


E= (7) 


and we see that the term A! in (6) is but a 
natural consequence of the hypothesis pro- 
posed in this paper. As to the factor 1/6 ez?, 
if we put e=5x 10", z=10~ into it, we get 


1 
6 BVA. 


1 
Gers x 102, 


which is again almost equal to the factor 
6x10? in the formula (6). 

Thus, so far as the above considerations are 
concerned, a consistent picture has been 
obtained about the relations among the earth- 
quake volume, aftershock area, and mechani- 
cal strength of the crustal material, and 
above all about the way how the mechanical 
energy is stored up within the earth’s crust 
until it is released from it in the form of 
seismic waves. 

According to the view presented in this 
paper, the energy of an earthquake is depend- 
ent mainly on the volume within which it 
has been stored up until the outbreak takes 
place. This represents a controversial view 
as compared with another theory that the 
energy of an earthquake is determined by the 
crustal strain x rather than by the earth- 
quake volume V, that is a large earthquake 
corresponds to higher stress and a smaller 
one to lower stress. 

This study was made while the writer was 
staying at the California Institute of Techno- 
logy, Pasadena, California. Valuable discus- 
sions on this problem which he could have 
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with Drs. H. Beniorr, C. Hewitt Drx, B. 
GutenserG, F. Press and C. Ricuter, all of 
the Institute, are gratefully acknowledged. 
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A Simple Method for Calculating the Correlation Coefficients. 


By 


Yoshibumi Tomopa 


Geophysical Institute, Faculty of Science, Tokyo University. 


Summary 


A simple approximate method for calculating the correlation coefficient 0 (x, y) between 


two probability variables w and y, say, is given in the present paper. 


The essential point 


of this method is to divide the whole range of each variable into two classes and denote 
the variable x(or y) by +1 or —1 according as it is larger or smaller than its mean 


value. 
-prove the relation (3) in section § 1. 


If the correlation coefficient between these new series of a and y is 7(x, y), we can 
By using this relation, we can calculate the required 


coefficient 9 (a, 7) with much less numerical labour than otherwise. 


Sak 


The correlation coefficient po (2, y) between 
two probability variables x and y, say, is given 
by 


poof +oo 
oe v=_ |_| Pe idedy, 0) 
Oxr0y —oco J —co 


where o,” and o,? are the variances of « and 
y respectively, and P(x, y)dzdy represents 
the probability that 2 and y are in the ranges 
x~(a+dxz) and y~(y+dy). To calculate the 
correlation coefficient according to this defini- 
tion, a rather laborious work is needed. 
Thus, it is usually considered impracticable, 
if not impossible, to calculate manually long 
series of correlation coefficients such as those 
needed for auto-correlation function. Here 
we are going to propose a simple method to 
calculate the coefficient. When the sample 
distribution is normal and the sample size is 
_large enough, the new method enables us to 
obtain a very good approximation of (1) with 
very little labour. The essential point of the 
present method is to divide the whole range 
of each variable into two classes only. The 
variable 2 or y which is larger (or smaller) 
than its mean value will be denoted by +1 
(or —1). Thus,.we get a new series of 
variables, all expressed by +1 and —1. The 
correlation coefficient 7 (a, y), say, for this 
new series is given by 
NeSN- 


1 (@;9) SL? (2) 


where WN. (or N-1) is the number of cases 
where z and y are the same (or opposite) in 
sign. The correlation coefficient (2) is very 
easy to calculate. As will be proved in the 
next section, o(2, y) and r(a, y) are related 
as follows: 


p(w, y)=sin4 7 (@, 9). 3) 


By using the above relation, we can easily 
calculate the required true coefficient 0 (a, y). 


SP 

The proof of the relation (3) is as follows. 
If the sample distributions of x and y are 
both normal, by using the normalized variables 


a and Y=", the probability P(X, Y) 


that X and Y are between X~(X+dxX) and 
Y~(Y+dY) is given by 


Ne 


PY 
1 — {X2-20XY-+ Y°} /2(1 — 07) 
~ on — 0" 
(4) 
Then, 
HH, “DW EN. 


== Nes Vas 8h | lux. Y) dX dY 


1st md 


8rd 4th 
quadrant quadrant 


_ 4 f. 5 -/ste— al 
mart poten ey: | 
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o(a, y)=sin - re 


$3. 
In this section, several numerical examples of the application of the present method will 
be given. 


Example 1. 
x: Students’ mathematics test marks y: Students’ language test marks 
Sample size: 128 
Regn ee ee ee ee 
Mathematics 
x | | i | | | 
| 35~ | 40~ | 45~ | 50~ | 55~ | 60~ | 65~ | 70~ | 75~ | 80~ | 85~ | 90~ | 95~ 
wi i aes a eee es Se) es eee ns — 
65~ 1 Ia’ ea | | 
70~ 1 jp Sze Sees | } 4 
75~ Pa | Sy OR te al | | 
Language 80~ | 3 8 | fi 5 | 4 | 1 | 
85~ ) peo atc gi) / a ota puee 1 
| | | | 
90~ | | 1 20 hte 5 3 4 
05~ | | | A ip h't 2 2 
o(true): 0.84 r(+1, -1): 0.66 p=Sin__ 77 Uso 
Example 2. Example 3. 
x: Students’ first test marks x: Weight of body (in kg) 
y: Students’ second test marks y: Height of body (in cm) 
Sample size: 100 Sample size: 125 
| 
| First test | Weight of body (in kg) 
7 i  # [ ‘| | Ms is os es ee “at a | ou | [ 
| |65~ | 70~ | 75~ | 80~ | 85~ | 90~ | 95~ | |40~| 45~|50~|55~| 60~| 65~|70~|75~ 
\¥ | ea eee | Ly | | | | 
}65~ | 1) = |150~| 2 be | 
9 | | 
Pe 70~ <i i | = 155~ | Sy | 2 | 3 2 | a 
aes | 2) 3) &) Cy ay > |160~ | 612" kien Greed | 
Sees | 8] 7] 8) 4 | 3 |165~ | 1 | Woya5 } 8) 8 | a 
3 85~ 4| 4 6 4 9 | 1 | S | 170~ | | Sales 3 4 3 a 
90~ | Lt 34 Sa te & |175~| | 3 lle f 2 
95~ | 1) 2 ees | 180~ | tect 
——EEEeE—————————————————————— = i } | 
op (true): 0.58 o (true): 0.61 
r(+1, —1): 0.36 r(+1, -—1): 0.31 
. 
p=sin5r: 0.54 p=sinr: 0.47 
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Example 4. 


w: Life-times of fathers y: Life-times of their children 


Sample size: 1072 


eae ee ey 
Fathers 
| = omit i plea fe estes P| g Pl =r TG a er | i | 
* 59 | 60 | 61 62 | 63 64 65 | 6 | 67 | 68 | 69 | 70 | 7 728 eas 
ret i oe | ‘ae i i srr i 
| or, | ee 5a) es | 
6 ey ig ec om 
62 BEE 256 Feil 6 Hele 2H asp aaitecs: | Sel ye «et of 
63 25) .25| 2.95) 2.9512. | 4. | 5. 2.75 1.25 | .25| 25 | 
64 1. 1.5 | 3,753. | 4.2518. | 9.25) 3. | 1.25 1.5 | 75) 1.25 | 
a5 (2. |1. |1. | 2. | 3.25) 9.5 [13.5 |10.76| 7.5 5.5/3.5 | 2.5 | 
66 | 52. | 2.25 5.25 9.5 10. ie.75 7.5 16. | 5.252. | 2.5 | 1. | | 
67 1,5 | 1.5 | 4.75) 3.5 /13.75| 9.75126.5 [25.75/19.5 /12.5 |13.75) 3.25) 5. | 1. | 
68 1. | 2. | 7.510. |10,2524.25)31.5 [23.5 [29.5 [13.25) 8.5 9.5 2.25 | 
E 69 | 5.25] 5. 12.75 18.25 16. 24. 29. 21.5 10. | 3.5 2.25 e 
a 70 Be | 2.5 5.7 18.7511.7519.5 22.5 19.5 14.5 6.25] 3.5 | 1. se 
71 | 3.25, 5. | 8.750.759. 14.7520.7510.35,8. |5. | 1. | 1 
72 | —“<}as|3. | 1.2817. | 7-6 bo.75)11.2510. | 8.5 2.75.5 
73 ine fa) 752751) 7.94, 6.5: 6.5 |. 7.8 | 6.25 3.25 5] .5 
| 74 1.5 | 1.5 | 5.25] 2.25| 2.5 | 6.5 | 3.25] 3.25 2 
| 75 | | | 12 3.5) 7.5 1.75] 5 
bee | eal ee aS mete 2a oe 
a7 | | | | 1.25) .25 | | 15) | 
bers: | | | | | | | 25 | .75 
zal | | | | | | 25] -25| 


A ™ 
pacinuie)e O-oll y(+1, —1): 0.36 p=sin-or: 0.54 
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Example 5. 


The autocorrelation correlograms of seismogram shown in Fig. 1 are computed by this 


new method and are shown in Figs. 2~5. 


AA athena IN a pclae ann 


flr) £0 
1.0 10 h 
5 5 
0| 0 
' / min. 4 F jmin. sy 
Fig. 2. Auto-correlogram for P wave Fig. 3. Auto-correlogram for S wave 
(Part 1) (Part 2) 
LO f,(e) 
1.0 10 
el ia 
‘| Ses MEANS: be 
— jmin 1 , [min ————s 
Fig. 4. Auto-correlogram for surface wave Fig. 4. Auto-correlogram for coda part 
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Correlogram Analyses of Seismograms by Means of 


a Simple Automatic Computer. 


By 
Keiiti Axt 
Geophysical Institute, Faculty of Science, Tokyo University. 


Summary 


An automatic relay computer for calculating autocorrelation coefficients, based on the 
simplified method of its computation proposed by Y. TOMODA, has been built for the pur- 


pose of analysing seismograms. 


Seismograms of about one hundred near earthquakes 


recorded at the Central Meteorological Observatory in Tokyo and at the Mizusawa Latitude 


Observatory were used for the analyses. 


Each seismogram was divided into several por- 


tions, and within each portion, the dependence of the prodominant period of seismic waves 
upon the magnitude of the earthquakes and upon the crustal structure of the wave path 


were studied. 


§1. 
What are the main factors that determine 
the frequency spectrum of seismic waves? 
This is a very interesting question from the 
standpoint of pure and applied seismology. 
The spectrum may depend both on the crustal 
structure at the origin of the seismic waves 
and also that along the passage through which 
the waves are propagated, as well as on the 
mechanism of wave generation at the origin. 
Although many studies have been made about 
the change in wave period with the epicentral 
distance, the variation in period with the 
earthquake magnitude, and so on, there still 
remain many questions to be investigated. 

In the present paper, the writer has made 
correlogram analyses of actual seismograms. 
~The correlogram. of a time series is the 
Fourier transform of the periodogram or the 
spectrum density of the series, and involves 
the same amount of information about the 
series as that of the spectrum density. 
Moreover, in many cases we can infer the 
spectrum directly from the correlogram with- 
out troublesome Fourter transformations. To 
make the computation process simpler, an 
automatic computer is urgently needed. Our 
computer is that based on the simplified 
method of computation of correlation coeffici- 
ent proposed by Y. Tomopa (1956). The 


Introduction 


computor is very easy to operate and main- 
tain. The cost of its construction is low 
and the high speed computation is possible. 


§2. The Automatic Computer and the 


Method of Analysis 
Let y() be a given time series, of which 
the mean value y(@) is 0 and the variance 


y*(2) equals to 1, and let the values of y(¢) at 
OOS Vie —— Ab AI im spoon pes ocO. Bia5, 
Y-1, Yo. Y1, Y2,°*:. Then the autocorrelation 
coefficient is defined by 
+7 
BAe = ae on 7 = Heel nl ( : ) 

A computer to calculate the autocorrelation 
coefficient according to the formula (1) needs 
the following parts; input devices such as 
tapes or cards, arithmetic elements to make 
additions and multiplications, and output de- 
vices such as a printing machine or some 
other indicators. So far as we stick to this 
regular method of computation, the comput- 
ing machine will become a large one, which 
cannot be very inexpensive. 

If we try to deduce the correlation coeffici- 
ent using the simplified method proposed by 
Y. Tomopa, however, the computer will be- 
come a very simple one. According to this 
method, the autocorrelation coefficient of y(Z) 
is given as 
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Fig. 1. 


Autocorrelation computer. 


Reet | 5 rhe)! 


where 7(k4¢) means 
(RAt)= lim —+- $* 
Y = lim Ss 
: ) nm > DH +1 r=a—% 
and sign(y,) is a function which takes the 
value +1 when the sign of y, is positive, and 
takes —1 when the sign of y, is negative. 
r(kAt) may be rewritten as 


N.(k)—N-(k) 
N..(k)+N_(k) 


where N.(k) is the number of sample 
for which y, and 7,4; have the same 
while N_(k) is the number of sample 
for which 7, and y+. have the opposite 
The number N of the total sample pairs 
equals to the sum of N.(k) and N_(k). On 
the other hand N,(0) equals to N by definition, 
so that if we get N(R) for k=0, 1, 2, --- , m, 
we have N-(k) by subtracting N.(k) from 
N,(0) and thus we can deduce the values of 
net) for k=O, 1 2. se gee by the use sor 


sign(yr) X SIgN(Yr+e) (3) 


1(kAt) = (4) 
pairs 
sign, 
pairs 
sign. 


the formula (4). 

Our computer gives automatically N+(%) for 
k=0, 1, 2, ---, 10 from the original record of 
y(¢). What we must know about y(t) are 
only its sign, so that the input and output 
systems and the arithmetic elements become 
very simple. 

In tracing the record to be analysed, the 
sign of y(¢) is read from it at regular time 
intervals, and when y(Z) is positive, a contact 
key is closed manually, while when y(Z) is 
negative, the key is made open. This con- 
tact key is connected to an input relay, de- 
signated by P, in the computer. Thus we 
can send input data into the computer direct- 
ly from the original record without any input 
devices such as a tape system. 

The multiplication in this computer is limit- 
ed to the following four cases: 


+1 Xl ek 
+1x-—l=+1 
=X alee 
—1x—l=+1 


and we can easily construct a circuit of re- 
lays for making this arithmetic. The addi- 
tion in this case is essentially a counting, and 
can be performed by means of an electro- 
mechanical counter. The results of counting 
of N.(k)’s are given at the indicator of the 
counter. 

To calculate N.(k)’s up to R=10 in one 
cycle of input data, a delay operation is 
necessary in our computer. Ten sets of two 
relays are provided for this purpose. The 
first two relays A, and B, are made to be in 
the same state at any moment as that of the 
input relay P, 4¢ before the moment. In: the 
similar way, the kth relays A, and By are 
made to be in the same state as that of P, 
kAt before the moment. 

Now we shall explain the details of our 
computer. A drum on which a seismogram 
to be analysed is wound as shown in Fig. 2 
is revolved with an appropriate speed. The 
reading of signs from the record of seismic 
waves is done through a cross wire, of which 
the intersection point is adjusted to coincide 
with the zero line of the recorded curve y(¢). 
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pulse time 
regulator. 


motor. 
Fig. 2: 


With the revolution of the drum, the record- 
ed curve y(¢) moves along the time axis, and 
when y(¢) is positive the contact of the key is 
closed manually. While y(¢) is positive, the 
contact is kept closed, and at the time when 
the curve passes through the intersection point 
of cross wire from the positive side of the 
zero line to the negative, the contact is made 
to open and kept opened as long as y(¢) is 
negative. The input relay P in the computer 
is connected to the contact of the key, and 
works just in the same way as the manually 
operated key. 

For the control of the operation of addition, 
multiplication and delay, there is a control 
unit composed of four relays. The four relays 


1,1, 2 and 2 of the unit work according to 
the time table shown Fig. 3 (left), in which 


Be lcat ae 
ae a oy a | 
iy ra ae 


Fig. 3. Time table of operation. 
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electromechanical 
relays. 


electromechanical 
counters. 


Ree~ DOWer source. 
b 


~input key. 


Se 


<i : 
seismogram. — crosg-wire, 


Computer with appurtenances. 


the line represents the duration of closed con- 
tact of each relay designated at the top, and 
the blank represents the duration of open 
contact. The unit time 4¢ can be adjusted 
from outside of the computer. We can 
obtain by combining the above relays a set 
of contacts of which the time table of opera- 
tion is as shown in Fig. 3 (right). In this 


figure, ~1 represents a contact acting re- 


versely to 1, i.e. the contact ~ 1 opens when 
the coil of the relay is energized, and closes 
when the coil is deenergized. For instance, 


the contact —1—~1— closes only when the 


contact 1 closes and 1 opens. 


p ~P 
Qj}; ~@) 6b; ~b1 


T 


~T 


Fig. 5. Multiplication and 


addition circuit. 


Fig. 4. Delay circuit. 
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A unit of delay circuit is shown in Fig. 4. 
When the contact of the input relay P is 
closed (opened) at the time of closure of 
—l1—~1-, the coil of the relay Aj is energiz- 
ed (deenergized) and the contact a; is closed 


(opened). While -i|- closes, the contact 


a, is kept closed (opened), thus the signal at 
the relay P is sent to the relay A; and delay- 
ed by 4¢. On the other hand, at the time of 


closure of —2—~ 2— the signal at the relay 
P is sent to the relay B,, and also delayed by 
4t. Further the signal at the relay A, is sent 
to the relay B, at the time of closure of 


—2—~2— and delayed by 4¢, thus the signal 
at B, at any moment is the one which was 
at the relay P, 24¢ before the moment. In the 
similar way, the delay of signal up to 104t 
can be performed by a connection of 10 such 
units. 

Multiplications and additions are done by a 
circuit of contacts as shown in Fig. 5. At the 


time of closure of —2—~2-—, the signal at 
P is multiplied by the signal at A;, and at 


the time of close of —1—~1-— the signal at 
P is multiplied by that at B,. An electro- 
mechanical counter C; counts the number of 
cases in which the signal at the relay Pis 
the same as that at the relays A, and B,, 
thus the indicator of the counter gives N,(1) 
in the formula (4). By the use of 9 other 
similar circuits of contacts and electromechani- 
cal counters, we can obtain N,(k)’s for k=2, 
3), Al ooo AKOy 

In the present analysis, the unit time dé 
has been taken as 0.5sec, so that from the 
data of sample size N, N.(k)’s for k=1, 2, 3, 
--+ , 10 are obtained within N/2 seconds. 

This computer was constructed by the Fuji 
Tsushinki Seizo K.K. (Fuji Communication 
Apparatus Mfg. Co., Ltd.), and was complet- 
ed in January, 1955. 


§3. Seismogram 


The seismograms used in this study are 
those recorded at the Central Meteorological 
Observatory in Tokyo and at the Latitude 
Observatory at Mizusawa, sixty from the 


former and fifty from the latter. The 
epicentres of all these earthquakes and the 
depths of focus for some of them are deter- 
mined by the Central Meteorological Observa- 
tory. Their magnitudes were independently 
estimated by C. Tsusor and H. Kawasumi 
(1952). 

The geographic distribution of the earth- 
quake epicentres recorded at the Central Me- 
teorological Observatory is shown in Fig. 6, 
and of those recorded at the Latitude Observa- 
tory is shown in Fig. 7. As shown in these 
maps, these earthquakes are grouped accord- 
ing to their epicenters into seven, in order to 
study the variations of statistical properties 
of seismograms with respect to epicentral 
areas from which they come. 


1») 
ox 
a 
‘ iy 
a G } 
eS x 
Fig. 6. Geographic distribution of epicentres of 


earthquakes recorded at the Central Meteorological 
Observatory in Tokyo. 


The Mizusawa Latitude Observatory has a 
history of fifty years of seimological observa- 
tions and there are abundant available data. 
Furthermore, this observatory is located just 
above the line of the very steep increase of 
gravity anomaly as revealed by C. Tsusor1 
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and his collaborators (1955). 
west of this line (named as the Morioka- 
Shirakawa line) is characterized by thick 
Tertiary formations which are developed only 
poorly in the eastern area. As the statistical 
properties of seismograms will depend on the 
crustal structure at the wave origin and also 
that along the wave path, it is expected that 
the correlograms of seismograms for the east- 
ern origin earthquakes may differ apprecia- 
bly from those of the western origin ones. 
The seismometer at the Central Meteorologi- 
cal Observatory is of the Wichert type with 


The area to the 


Fig. 7. Geographic distribution of 
epicentres of earthquakes recorded at 
the Latitude Observatory at Mizusawa. 


the static magnification 80 to 100, the free 
period 3.5 to 5.5 seconds and the damping 
ratio 2.8 to 5.7. Because of the magnification 
characteristics of the seismometer and the re- 
cording paper speed, the range of reliable 
period is estimated to be between 0.5-0.7 sec 
and 5.0-7.0 sec. 

The seismometer at the Mizusawa Latitude 
Observatory is the Omori type with the static 
magnification 20, the period being 31-38 sec. 
with no damper. 


§4. Results 

More than one hundred seismograms are 
loaned from the above two laboratories and 
put into the relay computer for correlogram 
analyses. We shall now enumerate some of 


the notable facts revealed by the analyses. 

1) The seismograms of earthquakes coming 
from the same or nearly the same epicentral 
area and with the similar magnitude have 
nearly the same correlograms when recorded 
at a station. 


Fig. 8. 


Correlograms of coda waves of earth- 
quakes occurring in the area PF. 


1.0f 


0.0r- SX 


a 


WS 


-1.0 


ig, 3o3 
quakes occurring in the area W. 


Correlograms of coda waves of earth- 


The correlograms of the coda parts of 
seismograms of earthquakes originated from 
the area EH and W are given in Fig. 8 and 9 
respectively. All of these earthquakes have 
the magnitude of 6.00.2. It is remarkable 
that each correlogram is of the similar and 
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rather simple form. By averaging the cor- 
relograms in each’ figure, we have a mean 
correlogram for each epicentral area. Both 
of the mean correlograms are of very simple 
types as shown in Fig. 10, so that-»we may 
assume the spectrum density of the following 


types. 

$(w) 

= we Bor +B?) 
(o—@,+28)(@—o,—78)(o+01+78)(o+a,—28) 


(5) 
Applying this formula to the two mean cor- 
relograms, we get 0.83zsec-! and 0.267 sec! 


1.0 


Fig. 10. Averaged correlograms, the thick line 
for the area # and the dotted line for the area W. 


c/s 


0.4 


Fig. 11. Spectrum densities corresponding to the 
averaged correlograms. 


0.8 


respectively as ; and § for the earthquakes 
originated from the area £, and 0.772 secu} 
and 0.19zsec-! for those from the area W 
respectively. The spectrum density is calculat- 
ed and is given in Fig. 11. 

In general, the correlogram of a part of a 
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seismogram does not always take such a 
simple form as above, but it is necessary for 
the reduction of data to assume a spectrum 
with a single predominant period for each 
part of a seismogram. We shall hereafter 
study mainly on the predominant period read 
from the spectrum or directly from the cor- 
relogram. 

2) How the predominant period changes 
with time after the beginning of P wave 
depends almost on the crustal structure at 
the epicentre and that along the wave path. 
In the case of a simple crustal structure, the 
predominant period is very short in the early 
part of seismograms and later increases to 
approach gradually to a steady value in the 
coda part. On the other hand in the case of 


100 200 


300 


sec 


Fig. 12. Change of the mean predominant period 
with time for earthquakes of the magnitude 6.0. 
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Fig. 13. Change of the mean predominant period 
with time for earthquakes of the magnitude 6.5. 


a complex crustal structure, the predominant 
period in the early part is nearly the same 
as that in the latter part. 

The change of average predominant period 
with time for the earthquakes of magnitude 
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Fig. 14. Change of the mean predominat period 
with time for earthquakes occured in the areas 
and W. 


Fig. 15. Correlograms of seismograms 
arranged in the order of the magnitude. 


6.0 and the epicentral areas B,, C and D, is 
shown in Fig. 12, and that for the earthquakes 
of magnitude 6.5 and the epicentral areas A, 


‘ 


Cc and D in Fig. 13. The earthquakes occurring 
in the area C are deep focus ones and the 
wave passage from there is regarded as of a 
simple oceanic type, while on the other hand 
the wave path from the area D is of a com- 
plex inland type, and that from the area A 
and B, is of an intermediate type. Fig. 14 
shows the period change in seismograms re- 
corded at Mizusawa for the earthquakes of 
magnitude 6.0 occurring in the areas E and 
W. As refered to before, the crust of the 
western area W of Mizusawa is more complex 
than that of the eastern area EF. The way 
of period change in the Figs. 12, 13 and 14 
agrees with the above statement. 

3) As to earthquakes having the same 
magnitude, the predominant period in the 
early part of seismogram shows a conspicuous 
variation with the place of epicentral area and 
wave passage, but the variation diminishes 
or dies out in the coda part as seen from the 
preceding figures. We may say, therefore, 
that the period of the coda wave is almost 
determined by the magnitude of earthquake’ 

4) As to earthquakes coming from the 
same epicentral area, the predominant period 
of every part of seismogram increases with 
the magnitude. But the rate of increase 
with the magnitude differs among parts of a 
seismogram, it is rapid in early parts and 
slow in later parts. 

The correlograms of the motions for 50 
seconds after the beginning of S wave are 
arranged in the order of magnitude in Fig. 
15. The largest magnitude is 6.8 and the 
smallest (the top one in the figure) is 5.2. 
All of these earthquakes occured in the area 
D. It is evident from the figure that the 
larger the magnitude, the longer becomes the 
period. 

To see the way of variation of the predomi- 
nant period with the magnitude, we plotted in 
Fig. 16 the period read from the correlograms 
of waves for 50 seconds after the beginning of 
S wave against the magnitude. The plotted 
points are scattered widely, but they may be 
divided into three groups according to their 
epicentral areas. The first group of earth- 
quakes is of inland origin such as those from 
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the area D, and shows the longest period. 
The second guoup is of the deep focus type 
with the shortest period. The last has the 
origin in the area A and B, and its period is 
intermediate. 


sec 


5.5 65 7.0 75 at 


Fig. 16. Variation of predominant period of S 
wave with the magnitude. 


6.0 


Assuming the following formula for each 
group 


logi9 T=const. + aM, (6) 


we have @ of about 0.3. The lines in Fig. 
16 are drawn for a=0.3. 

H. Honpa and H. Ito (1939) found a relation 
between the period of P waves and the 
magnitude of earthquake, and K. Kanal, K. 
Osapa and S. YosuizaAwaA (1953) extended the 
study. They got 7?=A between the period 
T of P waves and the amplitude A reduced 
to the place 100km from the origin. The 
magnitude M can be written as M=log,A+ 
const., so that we can convert the above rela- 
tion to the form as (6) and obtain w=0.5. 

After Kanat, the relation between the 
period of S wave and the reduced amplitude 
A is the same as that for P wave, but in 
the case of S wave his diagram seems to 
suggest the value of @ somewhat smaller 
than 0.5. 

According to the study of Hawaiian local 
earthquakes by A. E. Jones (1938), the upper 
limits of amplitudes of P and S waves of 
local shocks are found to vary with the cube 
of the period. In this case, the value of a 
becomes 1/3, and approximately agrees with 
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ours. 

As stated before, the predominant period 
of coda waves is mainly determined by the 
magnitude of earthquake. The periods of 
coda waves are plotted against the magnitudes 
of earthquakes in Fig. 17, in which earth- 
quakes of the epicentral area D are omitted. 
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Fig. 17. Variation of predominant period of 


coda wave with the magnitude. 


Fitting the same formula as (6) to this case, 
we obtain 


logio T= —0.82+0.22 M 


In Fig. 17, the periods of coda waves of 
local earthquakes in California are plotted 
against the magnitudes. These values are 
given in a paper by P. Bygerny (1947) who 
analysed numerous seismograms obtained by 
a seismometer of the Wood-Anderson type. 
The points are located on: an extrapolated 
part of our formula. 

B. GuTenperae and C. F. Ricuter (1942) 
obtained the relation logy 
between the periods T of the wave group of 
the maximum amplitude and the magnitude 
M. The value of the factor a is the same 
as in the case of coda waves, but the value 
of the constant in the formula is different. 
Thus the period of the coda wave is 4.8 times 


of that of the wave group of the maximum 
amplitude. 


$5. Acknowledgement 


The writer is very much obliged to Dr. W. 


M 


—1.5+0.22M 


Correlogram Analyses of Seismograms by Means of a Simple Automatic Computer 79 


Inouye, Mr. Y. Iwar and other officials of the 
Seismological Section of the Central Meteoro- 
logical Observatory and to Mr. T. Sucawa 
and other officials of the Mizusawa Latitude 
Observatory, who kindly lent him valuable 
seismograms recorded at the observatories or 
helped him to copy them. 

The writer’s thanks are also due to Mr. T. 
IkepaA and other staff members of the Fuji 
Tsushinki Seizo K. K. (Fuji Communication 
Apparatus -Mfg. Co. Ltd.) for their helpful 
advises given him about the design of the 
computer. 

The writer is greatly indepted to Prof. C. 
Tsupor for his guidance and encouragement. 

This study was made with the help of a 
Grant from the Department of Education. 


References 
AKAMATU, K.: 
1956 in this issue. 
BYERLY, P.: 
1947 The Periods of Local Earthquake Waves 
in Central California. Bull. Seis. Soc. 
Amer. 245.208. 


GUTENBERG, B. and RICHTER, C. F.: 
1942 Earthquake Magnitude, Intensity, Energy 
and Acceleration. Bull. Seis. Soc. Amer., 
2, loss 
HONDA, H. and ITo, H.: 
1939 On the period of the P waves and the 


magnitude of the earthquake. Geophys. 
Mag. 13, 155. 
JONES, A. E.: 


1938 Empirical Studies of some of the Seismic 
Phenomena of Hawaii. Bull. 
Amer., 28, 313. 

KANAI, K., OSADA, K. and YOSHIZAWA, S.: 

1953 The Relation between the Amplitude and 
the Period of Earthquake Motion. Bull. 
Earthq. Res. Inst., 31, 44. 

The Central Meteorological Observatory, Tokyo: 

1952 The Seismological Bulletin of C. M. O., 
Japan for the year 1950, 99. 

TsuBol, C., JITSUKAWA, A. and TAZImA, H.: 

1956 Gravity Survey along the lines of Precise 
Levels throughout Japan by means of a 
WORDEN Gravimeter, Part 7, Tohoku Dis- 


Seis. Soc. 


trict. Bull. Earthq. Res. Inst., Suppl. Vol. 
4, 311-406. 
TOMODA, Y.: 


1956 in this issue. 


JOURNAL OF PHYSICS OF THE EARTH, VOL, 4, No. 2, 1956 81 


TOMODA’s Method for Calculating the Correlation Coefficients 
as Applied to Microtremor Analysis. 


By 
Kei AKAMATU 


Geophysical Institute, Faculty of Science, Tokyo University, Tokyo. 


Summary 


ToMODA’s method for calculating the correlation coefficients has been applied to micro- 
tremor analysis and it is shown that the method can safely be used for the purpose, if 
the sample size is larger than about 300. 


A very simple method for calculating the 
correlation coefficient between two variables 
x and y has been proposed by Y. Tomopa 
(1956). The method consists essentially in 
comparing the algebraic signs of 4z and the 
corresponding 4y, where 4x and dy are the 
deviations of 2 and y from their respective 
mean values. If the number of pairs of x 
and y for which 4z and dy are of the same 
sign is N. and that for which they are of the 
opposite sign is N_, the correlation coefficient 
op between x and y is given by 


inane 2) 
p=sin|(- sae \\ 


correlation coefficients have to be calculated 
as in the construction of an autocorrelation 
correlogram. 

In the present note, this method of calcula- 
tion will be applied for microtremor analysis. 
By microtremor is meant here the so-called 
ground noise or ground unrest which have 
periods less than a few seconds. The seismo- 
grams used for the analysis are those obtained 
at Tokyo by two seismometers, one a 1c/s 
horizontal and the other a 3c/s vertical. 

Examples of the records used are shown in 
Fig. 1. On these records, time intervals of 
6-8 seconds were chosen, and the mean value 
within the interval was established. The de- 
flection of the curve /f(¢) for every 1/50 of 
a second (t)) was read as either +1 or —1l 


NEN 2 


This method is conveniently used especially 
in those cases in which a large number of 


BWV VV [Worn 


ee 
| Sec, 


Fig. 1. A 3c/s vertical seismometer B 1e/s horizontal seismometer 


according as it is above or below the zero 
line. Thus from 300 to 400 readings of +1 
were obtained. From these readings, 


1 noe gt! 
=—= AA 
b(t) V62V 62N a3 vay , 


has been counted, for various values of t=kro, 
k being.taken to be 1, 2, 3, ---, 30 for the 1 c/s 


seismometer and 1, 2, 3, ---, 20 for the 3c/s 
one. o,2 and a,” are the variances of 4a and 4y. 

The autocorrelation correlograms, which are 
the curves of #(r) plotted against t, are shown 
by broken lines in Fig. 2. The solid curves 
shown in the same figure were deduced by 
the customary method of calculation. The 
latter calculations would have been extremely 
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Fig. 2. Correlograms, 1c/s horizontal seismometer 


tedious if done manually. Actually, a relay were actually deduced according to 
computor FACOM belonging to the Fuji Com- 
munication Apparatus Mfg. Co. could be used 6™ 
for the purpose through the courtesy of that 
company for which the author expresses her 
thanks. (AKkamatu: 1956). 

The curves of solid and broken lines, agree 
well in general tendancy, although in some 
cases the correlation coefficients differ notably, 
even at the first significant figure. 

Based on the above results, power spectra 
of the microtremors were futher constructed. 
In order to get rid of the effect due to minor 


fluctuations in ¢(t) as seen in Fig. 2. e717! 2 10 < 
was multiplied to the original ¢(c), a being FYs 
tentatively taken to be 6-8. Thus the spectra Fig. 3. a. Spectrum from Fig. 2. a. 


TOMODA’s Method for 
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Fig. 3. b. Stpectrum from Fig. 2. b. 
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Fig. 3. d. Spectrum from Fig. 2. d. 
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FY 
Fig. 3. c. Spectrum from Fig. 2. c. 
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Fig. 3. e. Spectrum from Fig. 2. e. 
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The equipment can be used for making enormous and complicated calculations very quickly 
through operation of electromagnetic relays. 
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The main components are: Easy for programming 
Manipulating key board for tape perforators High speed in calculation 
Control board Very accurate in operation 
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Studies on the Instability of the Layer of Fluid Heated from Below. 


By 


Yoshinari NakaGAwA 
The Enrico Fermi Institute for Nuclear Studies, The University of Chicago. 


Summary 


The details of the experimental studies on the instability of a layer of electrically con- 
ducting fluid heated from below subject to an external magnetic field are reported. A 
summarizing discussion on the theoretical studies on the instability of the layer of auras 
heated from below under various circumstances is also presented. 
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§1. Introduction 


The thermal instability of the layer of vis- 
cous fluid heated from below has become the 
subject of theoretical investigation since the 
first attempt of Lord RaytetcH (1916) on the 
experimentally observed cellular convection 
under such circumstances by H. Benarp (1900). 
RayLeicH treated this problem by the method 
of perturbation, and showed that cellular con- 
vection of specified dimension should ensue 
when a non-dimensional parameter called the 
RAYLEIGH number exceeds a certain determinate 
critical value. As RaAyLeicu’s theoretical work 
assumed only boundary conditions of both free 
plane bounding surfaces of the layer, H. Jrr- 
PREYS (1926, 1928, 1930) has extended the theo- 
retical treatment to other bounding surfaces 
such as two rigid planes, one free and one 
rigid planes. Thus, Jerrreys has obtained three 
different critical values of the RAYLEIGH number 
under each set of the boundary conditions; 
651 for both free surfaces, 1108 for one free 


and one rigid surfaces, and 1706 for both rigid 
surfaces. 

A. R. Low (1929) discussed this problem 
from the standpoint of the energy balance in 
the fluid, and showed that the similarity of 
the theoretical treatment of this problem to 
that of the dynamic instability of a viscous 
flow between concentric rotating cylinders. 
The latter problem was theoretically investi- 
gated by G. I. Taytor (1923), and he obtained 
the criterion of instability in the similar: criti- 
cal value of a non-dimensional number. The 
critical values of the two criteria closely agreed 
with each other and Jerrreys obtained 1706 
under both rigid boundary conditions while 
Taytor obtained 1709 for his problem. 

A. PELLEW and R. V. SouTHWELL (1940) have 
summarized these earlier theoretical studies 
with discussions on the validity of exchange 
stability, boundary conditions and the charac- 
teristic number @?. Then, they introduced the 
method of variation in the evaluation of the 
critical value of the RaAyLeIcH number. Later, 
this method has been developed by S. CHANDRA- 
SEKHAR into a more refined form, and he has 
studied the instability of the layer of fluid 
heated from below under various circumstances 
(1952a, 1952b, 1953b, 1953c, 1954a, 1954b, 1954c). 
He has also applied the same method of vari- 
ation to other problems of the instabilities in 
hydrodynamics and hydromagnetics (1953a, 
1953d, 1954d). Thus, the theoretical treatment 
of this problem practically leads the general 
theoretical studies of the problem of instability. 
While CHANDRASEKHAR has been developing 
such theoretical treatments, independently W. 
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B. Tuompson (1951) studied the convection 
under a magnetic field, and the author (OG 
Naxkacawa and P. Frenzen 1955) studied the 
cellular convection of rotating fluids. How- 
ever, these theoretical studies were completely 
covered by the works of CHANDRASEKHAR (1952a, 
1953b). 

The results of earlier theoretical studies, 
especially Jerrreys’ results, have been examin- 
ed experimentally by R. J. Scumipt and S. W. 
Minverton (1935), R. J. Scummt and O. A. 
Sanpers (1938) and others. They confirmed 
that instability really sets in as cellular con- 
vection of a specified dimension when the 
critical value of the RAYLEIGH number is ex- 
ceeded. However, some disagreements between 
experiment and theory were noticed in the 
experiments studied by K. CuHanpra (1938), 
and QO. G. Sutton (1950) gave a theoretical 
discussion on this aspect. Sutton explained 
that these disagreements were caused by the 
presence of the other mode of convection 
which should be able to appear in the layer 
when the depth of the layer is less than the 
thickness of the boundary layer of fluids. 
Then he successfully interpreted the results 
of CHANDRA’s experiments. 

The results of the theoretical studies on 
the rotating fluids were confirmed experi- 
mentally by Y. Naxacawa and P. FRENZEN 
(1955) with respect to the theoretical criteria 
of the exchange stability. The other aspect 
of the theoretical result, that is, the onset of 
oscillatory convection, was firstly demonstrated 
by D. Fuitz, Y. Naxacawa and P. FRENZEN 
(1954) in a trial experiment with a layer of 
rotating mercury. Later, experiments on this 
aspect performed by the author (D. Futrz and 
Y. Nakacawa 1955) obtained a good agreement 
with the theoretical criteria especially prepar- 


ed for these experiments by S. CHANDRASEKHAR 
and D. D. Etpert (1955), 


Experiments on the instability of the layer 
of electrically conducting fluids under external 
magnetic fields have been performed by the 
author (Y. Nakagawa 1955)* and K. Jrruow 
(1956) and they have provided confirmations 
for the results of theory. In this paper, the 
details of the author’s experiments are reported 


with a general discussion covering the whole 
aspect of the theoretical studies of this problem 
under various circumstances. 


$2. The Equation of the Problem 


The equations of this problem are discussed 
by means of conventional vector operations*™. 
By including the influences of rotation and the 
magnetic field, we have the equation of motion 
in the following form which is called the 


equation of magneto-hydrodynamics: 


pv + o(V-7) V+ 2Qx V)—JX wl 


=—pP+p = y(V-V)+o»p?V+opU, (1) 


where V denotes the velocity, o the density, 
mw the magnetic permeability, P the pressure, 
vy the kinematic viscosity and U the gravita- 
tional potential. The term 2(2x V) denotes 
the Cortouis’ forces*** and the term Jx nH 
denotes the mechanical force exerted by the 
magnetic field H, to the electrically conducting 
fluids carrying the current density J****. The 
term due to the space charge is neglected in 
the equation because the motion of the fluid 
is naturally much smaller than the velocity of 


The equation of continuity becomes in the 
form: 
Op 

-(0V)=0:. 

9 au (eV) 

The electro-magnetic phenomena of the 
fluid are described by the MaxweELu’s equation 
of the following forms: 


ona OH. 


(2) 


* A brief note was published in Nature. 
** See S. LUNDQUIST (1952); some operational 
relationships are listed. 
*“* The justification from microscopic view 
points is given by the author in the foregoing paper. 
week See S. CHAPMAN and T.G. COWLING (1939); 
the interpretation of this term from microscopic 
view point can be easily obtained. 
ek Discussions are given by G.K. BATCHELOR 


(1950), S. LUNDQUIST (1952), E. BULLARD (1955) and 
others. 
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VX H=4z7J , (4) 
Vy-H=0, Gay 
y-E=(, (6) 


J=o(E+VX uA) , ye 


where E denotes the electric field, o the elec- 
tric conductivity, and the last equation is the 
Oxw’s law for the moving fluid with velocity 
V. The electro-magnetic cgs. units are used 
with the approximation of slow motion of 
fluid compared to the velocity of light. The 
thermal and kinematic phenomena of the pro- 
biem are related by the equation of energy”: 

Vie P ; F ) 

at =ECy Pete [= wet al A 
where T is the temperature, J the mechanical 
equivalent of heat, C, the specific heat for 
constant volume, « the thermometric conduc- 
tivity, and @ is the “dissipation function ”. 
Finally, we have the equation of thermal 
expansion, 


(8) 


o=p,(1—a@4T) , ey) 


where a denotes the coefficient of volume ex- 
pansion and 47 is the difference between the 
temperature at which the density is » and the 
temperature at which the density is gp. 
Before the onset of convection, with the 
heating from below, a state of thermal equili- 
brium should be attained in the layer of fluid 
due to the conductive heat transfer. This 
leads to the assumption of the initial state of 
being steady and relatively at rest. Then we 
obtain the following relations as the initial 
conditions: " 


yP=opU, (10) 
iO). (11) 
J=0, (12) 
E=0, (13) 


where the equation (10) denotes the hydrosta- 
tic relation and the equations (12) and (13) are 
ybtained with an additional assumption of 
iniform external magnetic field H. Integrat- 
ng the equation (11) with the assumption of 
miform temperature in horizontal direction 
ind higher temperature at bottom, we have 


; oT 
Totgeae. ( p=", ) (14) 
where 7, denotes the temperature at the bot- 
tom of the layer and # is a constant and 
denotes the initial vertical adverse temperature 
gradient. 

Now for the investgiation of instability, we 
assume that the quantities concerned with the 
unstable motion, i.e. the convective motion to 
be small enough so that squares and products 
of such quantities are negligible. Then, taking 
perturbation of the original equations (1) to 
(8) we have the following equations: 


oy 02(2 xX u)—j X vo 


=—pb+0? PP-u)+orruty'PU, (5) 


00’ 
py FU P)e+(P-u)p=0 , (16) 
___ Ok 

Vxe= Maes (17) 
yxh=4rj, (18) 
V-h=0, (19) 
V-e=—yH- (7 Xu), (20) 
j=o(e+ux A) , (21) 

00 P 
—(u-k)B+ - -u)=Kp70 22 
ae (u-k)B+ alc, ? u)=«Kp6 , (22) 


where k is a unit vector in the vertical direc- 
tion and u, j, Dp, 0’, e, h and @ denote the 
quantities concerned with the convective 
motion. 

Rewriting the equation of thermal expansion 
(9) by these quantities, we have 


0 =—eao , (23) 


“where o’ is negative as the density decreases 


with the increase of the temperature. The 
equation (23) suggests that due to the factor 
a the magnitude of the density variation be 
always smaller than the magnitude of the 
temperature variation**. Therefore, to the 


* See S. GOLDSTEIN (1943) p. 606. For gas Cy 
can be replaced by C, with « replaced by « Cy/Cp. 
(see G. K. BATCHELOR 1954). 

** This reason also justifies the elimination of 
terms in the perturbation equations which involved 
the variation of p caused by the initial temperature 
gradient in the fluid. See the discussion given by 
A. PELLEW and R. V. SOUTHWELL (1940). 
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same order of approximation, we can simply 
write the equation (17) as 
y-u=0. (24) 
Thus, using the equation (24), the equations 
(15) and (22) become 


pots 02(2xu)—j X ui 
=—pp+prpu—paopu , (25) 
oF (ub) = 9°0 (26)* 


Eliminating e in the equation (17) by equations 
(18) and (21), we have 


Oh 


——(H-p)u=np7h , (27) 
aL (H-7)u=77 
where 7= _ . Similar elimination of j in 
4nou 
the equation (25) by (18) gives 
Ou Dp 


+ 2(2 x u) ep hh x H=—y— 
4ro 0 
+yp?utrok , 


Ot 
(28) 


where ry=ga@. Operating curl (yx) to the 
equations (27) and (28), we obtain 


0 » 
OP (HDB 


(29) 
OA 


A +2(2-P)u—7 


Ip ED B= A417 x b)O 
(30) 
where A=pxXu and B=pxh. Another opera- 


tion of curl to the equation (30) gives 


© (py) +.2(0-p)A —"-UE- p)prh 
Aro 


Ot 
=vyp'u—rpxVxk)é , (31) 
where the following relations are used 
VxA=—pu | (32) 
and yxB=-frh. } 


The equations (26), (27), (29), (30) and (31) 
give a complete set of equations for further 
investigation of this problem, but for the 
simplification and convenience, we will confine 
the discussion to the vertical components of 
the various vector quantities, u, h, A and B. 
Then from the above set of equtions we have 


a ae 
(Gr \p= Bu . (33) 


& —ap" \=Fu, (34) 
(G1 )O=FE (35) 
(G7 k= Bw oases 39) 


| 
(G7 — )w= EC + ge rent rr? 1) 
where w, %, € and @ denote the vertical com- 
ponents of u, h, A and B respectively, and EZ, 
F and p.2 are the scalar operators defined as: 


E=2(8-7), + 
F=(H-9), (38): 
ee o2 

and Vi = ae tay : 


Eliminating x between the equations (34) and. 
(37), we have 


0 2 gi 2 es 2 2 
ee * Ne 48 yee ir ie 


@- bl Pecos 2 20 9), 
+E(5, w= ro mr?) . (39) 


Similarly, eliminating ¢@ between the equations. 
(35) and (36), we obtain 


(Cron) ek 


=E(5,—1* wo , 


Finally, eliminating € and @ in the equation 
(39) by using equations (33) and (40), we have’ 
an equation for w alone: 


Bos es\ Sf OV Oke SE Ven eee 
& rm )t(5, w\(a 10) — gag pre 
6) (6) : 
Ee — 2 ae 72 
: Gs (5, 9?) w 


(Ww MG) 


<P rw . (4) 
Azo 


(40)! 


The equation (41) is the basic equation for 
the investigation of the stability and the same 


* The similar equation results for compressible 


fluids: see the studies of H. JEFFREYS (1930) and Y.. 
NAKAGAWA and P. FRENZEN (1955). 
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equation is also satisfied by @, ¥, € and ¢ as 
it is easily seen from equations (30), (31) and 
(32). 

Under other circumstances, such as when 
no rotation or no magnetic field is presented 
in the problem, a simpler form of equation 
resulted as the basic equations which can be 
easily derived from the equation (41) by elimi- 
nating unnecessary terms. 

§3. The Functional Form of Solutions and 
the Boundary Conditions 


»The experimental observations of cellular 
convection approved the first step of the 
separation of variables z and y from z in the 
solution. Also, the pattern of symmetry of 
cells suggested a further speculation on the 
function x and y. A discussion was given by 
RAYLEIGH on the analogy of the vibration of a 
horizontal membrane, we can write as f« p,2f, 
where f denotes a function of « and y. The 
customary expression of the above relation in 
this problem is as follows, with the solution 
of the form of f(z, y)g(z)I(t), 

dpff+af=0, (42) 
where d is the depth of the layer and a? is 
the “characteristic number”. The important 
meaning of the “characteristic number”, @’, 
has been discussed by various authors (see 
Introduction). The practical meaning of @ 
is easily seen in its relation to the size of the 
cell. For instance, consider a hexagonal cell 
of the side length Z, then f is given by the 
following function discovered by S. Curisto- 
PHESON (1940): 


» fa} 00s ae 3at+y)+ bas ar Lh 3 «—y) 


4nx | 

—- ; (43 
+ Cos 37,4 f ) 
where f,) is a constant and » is an integer. 
Substituting (43) in (42) we have 


(44) 


Consequently, @ represents the horizontal 
wave numbers of the cellular convection. 
Similar relations as (44) are easily obtained for 
other specified shapes of cells such as squares 


or rectangulars. 
The boundary conditions which are immedi- 
ately apparent to be satisfied by w are 


w=(0) (45) 


at the top and bottom surfaces of the layer. 


Then assuming the constant temperature of 
boundary surfaces, we have 
d= (46) 
at the top and bottom of the layer which can 
be justified under the usual circumstances*. 
Further boundary conditions are obtained 
from the horizontal motion of convection. As 
no horizontal motion becomes possible at the 
rigid surface, the following conditions result 
from the equation of continuity, 
Ow 
Oz 
on rigid surfaces. Similarly, from the maxi- 
mum horizontal motions at free surfaces the 
following conditions result 
Pw jee 
Oz? Oz 
on free surfaces. 
The boundary conditions for electromagnetic 
quantities are 


=0 and €=0 (47) 


=0 and (48) 


Cee — andy — 0) (49) 
on a perfectly conducting surface and 
e—0) (50) 


on a free surface adjoining a vacuum. The 
equations (49) and (50) are derived from the 
theory of electricity and magnetism and have 
been introduced by CHANRASEKHAR (1953a, 
1954c). 

The above conditions are the basic boundary 
conditions, and in order to define the unique- 
ness of the solution the same number of con- 
ditions as of the order of the basic equation 
are derived in the following manners. For 
instance, from equations (33) and (37) we can 
derive the following conditions by substituting 
ee) and aie bespecvely 


* The cases such as non-conductive tenet 
00/02=0 and others are excluded because they are not 
essential points of the study. See’G. K. BATCHELOR 
(1954), H. JEFFREYS (1926), A. PELLEW and R. V- 
SOUTHWELL and others. 
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a 
ep =O, 51 
G xv) ou 


(5,27? pr EE 4 Fy ~=0, (52) 
ot Aico 

where both conditions hold on the top and 
bottom boundary surfaces. In this manner, 
the necessary number of boundary conditions 
are easily obtained even for the cases when 
the high order of partial differential equation 
is involved in the study. 


$4. The Exchange Stability and the Over- 
stability 

The time dependant function /(¢) in the solu- 
tions is customary assumed as an exponential 
function, 

I(t)= er" , (53) 
where ~ could be a real or complex number. 

The principle of the exchange stability as- 
sumes the reality of p as discussed by Jrr- 
FREYS (1926) and the criterion of the instability 
is defined by the relation, 

p=. (54) 
While the principle of overstability (the name 
“overstability” is used after A. S. Eppincton 
1926) assumes the complex value of ~ so that 
the criterion of instability is defined as 

br=0 (55) 
where p=p,+7pi, and p, and ~; denote the 
real and imaginary parts respectively. 

In other words, the principle of exchange 
stability assumes simple monotonous growth 
or decay of perturbations and the criterion 
indicates a steady onset of perturbed motion. 
The principle of overstability assumes the 
oscillatory motion of perturbations and con- 
siders the growth or decay of the amplitude 
of such oscillations. Then, the criterion indi- 
cates a steady oscillation of perturbed motion. 

The validity of these two principles of the 
stability have been discussed by several 
authors. When the solution of the problem is 
known, the validity of each stability is discus- 
sed from the so-called frequency equation of 
the stability. Examples of such discussions 
have been given by CHANDRASEKHAR (1952a, 
1953b) and the author (Y. Nakacawa and P: 


FRENZEN 1955). Another type of discussion is 
also presented by CHANDRASEKHAR (1952a) with 
the reference to PELLEw and SouTHWELL. 

A summarizing discussion on the validity 
and the dependence of the RayLeIcH number 
upon the important non-dimensional numbers 
is presented below covering the whole aspect 
of the problem under various circumstances. 
Before entering this discussion, we will intro- 
duce the conventional notation D which was 
defined as 

0 0 a2 

D=—= d— C2 

0z Oz o Sam d 

Then, we will write the solutions in the fol- 
lowing forms: 

w= f (x, y)Wia)er’ , 

A= f(x, yOZ)e”’ , 

C=f (x, y)Z(Z)e”' , 

1=F (w, YX@)e , 

d=f (a, yY@e" , 

where d@’pif (x, y)+a@f(x, y)=0. Further, con- 

fining the discussion to the case of @ and H 

both in the vertical direction, we have the 

following expressions for the scalar operators 

E and F, 


(56) 
rf 


(57) 


E=2 
d D (58) 
and roe D (59) 


The substitution of 2 given by (57) into the 
equation (41) gives the following equation 
after a little calculation; 


(ann Poo 
{PO — (oa) |p (Daw 


+ 7{2® _(p»_ayl {P& _(p»_ ay" pw 
K u] 
pd? 


#400 
x ie ~(P—a)| ap |W, 60) 

where 
Ie! yo (61) 


Ov 
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is the magnetic number*, 
pat 


a (62) 
the TayLor number*, and 


Rees 
KY 


(63) 
is the RayLeicH number. If we write 


en 


_{Pa-(D'=@)}[(wa-(D*— a) q- (D°~ @)) - QD*(D= a) W + T{Pq- (D8 @?)){q- (D?= a) DPW 


y 
q=— , ’ 
Vv Kk 


(64) 


] 
oOo =— 
H 


then by combining the expressions, we have 


qd 2 
PEE altgg Les (65) 
K | 
where P is the PranptL number. Snbstituting 


these notations in the equation (60) and rear- 
ranging the equation, we have 


a{ oq (D?— a*)}[{wq- (D?- a?) }{q-(D?- @)}-@D|W 


In this equation (66), we can see the compli- 
cated dependence of the RayLeIcGH number 
upon the various non-dimensional numbers of 
the problem. Under the circumstances con- 
sidered, the RAYLEIGH number is a function of 
qa. @, Q, Tf and a® as well as the actual 
functional form of W. 

Now, dividing the problem into the follow- 
ing cases, we will discuss the details in each 
of the specified cases: 

(i) Q=0, T=0, (the case of RayLeicu, 
JEFFREYS and PELLEwW and SouTHWELL), 

(ii) @=0, T=0, (the case of fluid subject- 
ed to rotation), 

(iii) QO0, T=0, (the case of electrically 
conducting fluids subjected to an external 
magnetic field), 

(iv) Qs<0, Ts<0 (the case of electrically 


- conducting fluids subjected to the simultaneous 


action of magnetic field and rotation) 
in ye Q=0)°7 =0; 
Let @ and T be zero in (66) and eliminating 
the same operators we have 
a ogee) DoW 
Cw 


p= 
(67) 
In (67) if we apply the principle of overstabi- 


lity, we can write 


_{Pq—(D'—@) }{q—(D?—@)} (D—@) W+-T {Pq (De —@)} PW 


piles (66) 


WEF - (68) 
The substitution of (68) in (67) indicates an 
imaginary value of R. As such R must be 
excluded, we can conclude that under such 
circumstances the principle of overstability 
must not be feasible for the criterion of 
stability. Now, putting g=0 in (62), we find 
that R is the function of a@ and W 
alae 
ew 
Therefore, the criterion based on the exchange 
stability can be uniquely determined by the 
specification of a? and W. Thus the criterion 
has been obtained as the minimum value of 
R with the specifications of @ and W in 
theories and which has been confirmed by the 
experiments (see Introduction). 

Hereafter, in the paper the difference of W 
for different boundary conditions and the 
reality of W are understood. Also, the defini- 
tion of the criterion as the minimum value of 
R with the specifications of a? for the ex- 
change stability, and the minimum value of 
R with the specification of a and 7 for the 
overstability are understood. 

(ii) Q=0, T0. 

In the similar way of obtaining (67) from 

(66), we have 


R= (69) 


(70) 


pitts a*{q—(D'—a")}W 


* More general forms of these numbers are Q= 


Sea Oe d? and T= 


AQ2 cos2 
sella where # denotes 
a 


Ov 


the angle of inclination of H or 9 to the vertical (S. CHANDRASEKHAR 1953b, 1954b). Also »“@ some- 
times called the HARTMANN’s number (R. C. Lock 1955) after the first work of J. HARTMANN (1937). 
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In equation (70), we can see the RAYLEIGH 
number now depends on the parameters P, 
q,@ and T. Also the feasibility of the over- 
stability can be seen in (70) as the odd order 
of q appears both in the denominator and 
numerator. Now, by putting q=0, that is, 
the criterion based on the principle of the ex- 
change stability, we find that the PranpTL 
number disappears from the equation with gq. 
Therefore, we can conclude that the criterion 
of the exchange stability is free from the 
PRANDTL number of fluids and the criterion of 
the overstability to the contrary. 

The criterion of the exchange stability has 
been obtained as a function of JT, and the 
criterion of the overstability has been obtained 
as a function T and P. Regardless of the 


pe _(Pa-(D—a)(oa— 


a{woq— 


Again by putting q=0 in (73), we conclude 
that criterion of the exchange stability is free 
from the PranpTtL number and w of fluids, 
while the criterion of overstability is depend- 
ing upon such numbers. Though the feasibi- 
lity of overstability can easily be seen in the 
equation (73), the overstability only becomes 
valid under astrophysical conditions as shown 
by CHANDRASEKHAR (1952a). With the assigned 
external conditions such as P, @ and ow, the 
general characteristics of the criteria become 
R.x«Q (when Q->co) (74) 
for both principles of stabilities (see the section 
7 for the experimental confirmation). 
(iv) @=<0, T><0. 
Then the expression (66) remains for the 


discussion. Putting g=0, we have 
Ie 
_{(D—a?!— QD" D?— a) W+ TD) W 
a{(D*—@)—QD*}W 
(75) 


From (75), we can conclude that the criterion 
of the exchange stability is still free from the 
PRANDTL number and w of fluids. However, 
the complicated characteristics of the depend- 
ence of R on Q and T in (75) have been reveal- 
ed by CHANDRASEKHAR (1954¢c). The overstabi- 


(D—a)}{q- 
(D?—a?)}W 


difference due to boundary conditions, the 


general relations of such criteria become 
R.« T3/2 

for the exchange stability and 
Rf l (when TP) (72) 


for the overstability, where the subscript c is 
used to denote the critical value of the RAYLEIGH 
number and superscript * is used to denote the 
criterion based on the overstability. The ex- 
perimental confirmations of these criteria and 
the theoretical condition of the validity of the 
overstability in terms of the PRANpDTL number 
have been obtained (see Introduction). 
(1) F=0; @=0. 

After eliminating the term including T in 
equation (66) we have 


(when T-> °°) (71) 


(D?—a?)}— -QD*|(D?— a?) W (73) 


lity is feasible as we can see from (66) and 
the criterion of overstability should include all 
the parameters. Regardless of the details of 
both criteria, the general characteristics of R 
and Q become 


Rex Q (when Q—> 0) 
for assigned values of T or T, P and w.* 


(76) 


§5. The Method of Variation 


The method of variation has been introduc- 
ed to this problem by PELLEw and SouTHWELL 
and has been developed into more refined 
manners by CHANDRASEKHAR as mentioned in 
the Introduction. Except for the case of the 
boundary conditions of both free planes, this 
method reduces a great amount of the practi- 
cal calculation of the criterion of the instabi- 
lity. Fortunately for the case when the 
boundary surfaces are both free planes, the 
form 


W=cos 7z (77) 


(where Z is measured from the center of the 
layer) becomes the exact solution of W under 
all circumstances, and we can easily evaluate 
the criterion by the suse aUaS of (77) in the 


* 


The experimental study is is now N under. prepara- 
tion by the author. 
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equations given in the previous section. Also, 
as the general characteristics of the functional 
relations between the RAYLEIGH number and 
the non-dimensional numbers in the criterion 
mostly remain for all boundary conditions, the 
solution (77) is an important tool for the in- 
vestigation (see the discussions, S. CHANDRA- 
SEKHAR 1952b, 1953b, 1954c, Y. Nakacawa and 
P. FRENZEN 1955). 

The practical criterion of this problem is 
the minimum value of the RayLeicH number. 
Therefore, if we can evaluate such minimum 
value of the RAYLEIGH number with the neces- 
sary specifications of @* or a? and 7, under 
certain assigned external conditions such as 
@, P, T and w, we need not have to solve the 
complicated equations of this problem. In 
short, the principle of the method of variation 
is based on such an idea and enables us to 
determine the criterion without obtaining the 
exact solutions of the problem. 

The detail of this method is illustrated below 
in relation to the following experiments for 
the case of Q*0 and 70. Under these cir- 
cumstances, as the overstability is only valid 
for astronomical conditions we will discuss the 
criterion of the exchange stability. From (73) 
by putting g=0, we have 


{(D—@)—QD*}(D?—a) W 
a*w ; 


Rewriting (78) we have the governing equa- 
tion of W in the following form: 


R=- (78) 


(2? —a)G+a0RW=0, (79) 
where 
G={(D—@yY—-QD}W. (80) 
The appropriate boundary conditions become 
W=G=0 on the boundary surfaces (81) 
and DW=0 ona rigid surface 
D?W=0 on a free surface , oe 


where the second relation of (81) was derived 
from (52) with the usage of (39), and this 
derivation can be easily justified. 

After multiplying by G and integrating (78) 
over the whole depth of the layer, we have 


-1/ 


+1/2 i 
| G(D'—aGdz 
—— aa ; (83) 
«| GWdz 
=/ 3 
where the origin of Z is taken at the center 
of the layer. The integral in the numera- 
tor and the denominator are easily transform- 
ed into the following forms after using partial 
integrations and the boundary conditions: 
the numerator, 
+1/2 
{(DG)?+@G}dz , 


(84) 


\y * G(D'—a)Gdi= { 


—1/2 -1/2 


the denominator, 

fi "CW \i- * ((D?— a?) Wy2-+Q(DW)']d3 . 
=1/2 -1/2 (85) 

With these relations (84) and (85), the equation 

(83) becomes 


|" (DG) +a'@)}dz 
R= = pei Ee 
a |" (D'—a) W}+QDW)"Nas 


-1/2 


, (86) 


which denotes the positive reality of R. Now, 
consider an arbitrary variation OW in the 
relation (86), then the effect of such variation 
on R to the first order of dW can be ex- 
pressed as 


sr=—_(0l,—aRoh) , (87) 
at, 
where 
=\" ” (Da) Wy+Q(DW)'|dz, (88) 
-1/2 
an=2\"" {(DG)(D8G)+.@GdG} da , (89) 
-1/2 
ano)" [{ (D2 a?) W }{(.D?—a?)0W } 
3/2 
+Q(DW)(DeW)dz, (90) 
and 6G stands for 
{(D?—@)?—QD}3W . (91) 


The expressions for 6/, and d/, can be rewrit- 
ten by partial integration and boundary con- 
ditions 
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DOW=0 on a free surface, | (92) 
DdW=0 on a rigid surface. 
After some calculations, the equation (87) be- 
comes 
Q (+1/2 = 
eae aA SG{(D?—a)G+a2RW } a2 . (93) 
Ato J-1/2 


In the equation (93) it is easily seen that the 
condition 6R=0 is only satisfied for arbitrary 
variations. in W provided 


(D?—a)G+@RW=0, 


that is, if W satisfies the original governing 
equation (79). Therefore, it follows that the 
exact solution of W would satisfy 6R=0 and 
if we evaluate R according to (86) with such 
a solution, we will obtain the critical value of 
the RAYLEIGH number as the minimum value 
of R. In practical procedures, the above 
principle leads us to the followings: assume a 
function G including one or two parameter 
and satisfy the boundary condition (81). Then, 
we can determine a function for W as the 
solution of the equation (80) satisfying the 
boundary conditions (81) and (82). Substituting 
such G and W in (86) and minimizing the 
value of R with respect to the parameters and 
@ (or a? and 7 if involved), we will have a 
good approximate critical value of the RAYLEIGH 
number. As W of a combination of even 
functions of 2Z satisfies both rigid conditions, 
we can choose the function G in such a form. 
Also a combination of odd functions of W of 
(77) satisfies the one rigid and one free condi- 
tions, provided the depth of layer, d, be replac- 
ed by d/2, the function G can be chosen in 
such a form under these boundary conditions*. 

For the present problem the critical values 
of the RayLEeIGH number have been obtained 
by CHANDRASEKHAR (1952a) after carrying out 
the calculation following the above method of 
variation (see the resultant relation (74)). 
Under other circumstances, as referred previ- 
ously, CHANDRASEKHAR has obtained the critical 
values of the RayLeicH number with appropri- 
ate choices of the function G convenient to the 
problems (S, CHANDRASEKHAR 1952a, 1953b, 


NAKAGAWA 


1954b, S. CHANDRASEKHAR and D. D. ELBERT 
1955). The values of R, thus determined by 
the successive approximations of adjustments 
of the parameters in the function like G, and 
a and sometimes a? and x have provided suffi- 
cient accurate values because the experimental 
results mostly have been obtained in the 
neighbourhood of such values (see Y. NakaGa- 
wa and P. Frenzen 1955, D. Fuuttz and Y. 
Nakacawa 1955 and the experimental results 
in this paper). 


$6. The Method of Experiment 


The difference of the amount of heat trans- 
fer through the layer between conduction and 
convection has been used as the base of the 
experimental determination of the critical 
values of the RayLEIGH number. When values 
of the temperature gradient are plotted against 
the amount of heat transfer or the amount of 
heat applied at the bottom of the layer, the 
onset of convection is easily detected as a 
change of the relations between these quanti- 
ties in a diagram. 

Such a diagram was firstly made by Scumipt 
and Minverton (1935), and they found a sharp 
break in the relation between the temperature 
gradient and the amount of heat applied at the 
bottom of the layer. By defining the critical 
value of the adverse temperature gradient at 
the point of break, they obtained the critical 
value of the RayLeIGH number in close agree- 
ment with the theoretical value. Later, ex- 
periments by Scumipt and Sanpers (1938) con- 
firmed the onset of the theoretically specified 
cellular convection at that point of break by 
optical observations. 

Though in most of these experiments, the 
temperature gradient was evaluated from the 
temperature measured at the top and the bot- 
tom of the bounding surfaces of the layer, the 
author has developed a method of direct mea- 
surement of the adverse temperature gradient. 
The latter method, especially measuring the 


Studies of the form of solutions of W under 
various boundary conditions have been made by 
H. JEFFREYS (1926), A. R. Low, A. PELLEW and R. 
V. SOUTHWELL. 


San 
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adverse temperature gradient inside of the 
layer has provided various important informa- 
tions. For example, the existence and the 
period of oscillation of the overstable convec- 
tion have been confirmed by this method of 
measurement (D. Futtz and Y. NaKaGAwa 
1955). In the following sections, the details of 
such method of experiment are discussed with 
the results of experiments for the case of a 


layer of mercury subject to an external magne- 
tic field. 


$7. The Experiments and Results 


The experiments are designed to investigate 
the functional relation between the critical 
value of the RAYLEIGH number and the magne- 
tic number Q, theoretically predicted by 
CHANDRASEKHAR (1952a). Excluding the over- 
stability under the present experimental con- 
ditions the relation of (74) becomes 


R.x«Q. (Q>1,000) (74) 
The general assembly and the schematic 
arrangements of the experimental apparatus 


are illustrated in the Fig. 1 and 2. As a 
permanent magnet* was used and it only 


MAGNET | 
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* The magnet was originally constructed ‘as the deflection magnet for cosmic ray experim 
very kindly made available for this experiment by Prof. MARCEL 


University of Chicago. 


provided a fixed intensity of the field, different 
depths of layer are used in the experiments. 
The pole surfaces of the permanent magnet 
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Fo eeR 1954 HA 
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Fig. 1. The front view of experimental set up. 


The schematic diagram of the experimental arrangements. 


ent and was 
SCHEIN of the Dept. of Physics, 
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were 17.8cmx2lcm rectangulars and the gap 
between the pole surfaces was 21.6cm. The 
intensity of the vertical magnetic field varied 
from the minimum of 1,375 gauss at the 
center of the gap to maxima of 2,250 gauss 
at both pole surfaces. Around the center of 
the gap the surfaces of the same intensity of 
the field formed coaxial hyperboric surfaces. 

A Pyrex-glass cylinder (A) of inside dia- 
meter 14cm and height 8cm was used for the 
container of mercury. The leveling screws 
(B) served to adjust the position of the cylinder 
in each experimental determination of depth. 
At the bottom of the cylinder a non-inductively 
wound resistance heater (C) was placed beneath 
a plastic coated copper plate. A five-couple 
copper-constantan thermopile (D) immersed in 
the mercury was used for the direct measure- 
ment of the internal adverse temperature 
gradient. The junctions of this thermopile 
were fixed at the levels about 1/2cm from the 
top and the bottom surfaces of the layer 
respectively. Also at each of the levels the 
junctions were spaced approximately 1cm 
apart along a diameter of the cylinder. Thus, 
the thermal e.m.f. of this thermopile elimi- 
nated the effect of the boundary layer of fluid 
and the variation of the temperature due to 
the relative positions of junctions and cells. 

A single thermocouple was used to obtain the 
estimate of the temperature of mercury. The 
measuring junction (E) of this thermocouple 
was placed at half way of the layer and the 
reference junction was placed in a constant 
temperature bath. To remove heat at the top 
surface, a cooled nitrogen circulation system 
was used, which also served to keep a const- 
ant average temperature of mercury. This 
system further worked to provide a_ better 
experimental approximation to the theory on 
the establishment of the linear temperature 
gradient through the layer. 

With these arrangements, experiments were 
performed with four different depths, namely, 
3, 4,5, and 6cm layers of mercury. At each 
determination, the positions of the cylinder 
were adjusted and the levels of the central 
plane of the layers were kept to coincide with 
the central plane of the gap of the magnet. 


Though this procedure reduced the variation 
of the intensity of the magnetic field to the 
minimum state, about 10% residual variation 
of the field was observed around the cylinder. 
Various rates of heating were applied at the 
bottom of the layer by the heater (C) with D. 
C. The current was adjusted by a combination 
of sliding resistance and measured by a D. C. 
ammeter. Then the rate of heating was de- 
termined by the square of the heating current 
to a sufficient accuracy. At each fixed rate 
of heating, a continuous record of the adverse 
temperature gradient was obtained by a system 
of d.c. amplifier and milliammeter recorder 
connected to the thermopile (D). The accuracy 
of this system was secured to the order of 
+1yV, and the system reproduced the record 
within the difference of the thickness of the 
line for independent trials of same rates of 
heating. The rates of circulation of cooled 
nitrogen were also adjusted in each experi- 
ment considering the rates of heating applied 
at the bottom of the layer. Thus, a number 
of records of the adverse temperature gradient 
was obtained for each of the four depths of 
the layer. These records were classified into 
typical three different types illustrated in Fig. 3. 

The records of type (a) were classified as 
the records of conductive heat transfer because 
the trend of the adverse temperature gradient 
were easily interpreted as the gradual esta- 
blishment of a stationary thermal equilibrium 
in the layer. This classification was support- 
ed from other standpoints. Thus, records of 
type (a) were only obtained in the lower rates 
of heating and no convective motion was 
observed when such records were obtained. 
The plotting of the final value of adverse 
temperature gradient (8;) against the corre- 
sponding rates of heating (J,) in a diagram 
gives the linear relation dd shown in the Fig. 
4. Considering the general relation of con- 
ductive heat transfer; 


K=RB , (94) 


where K denotes the amount of heat transfer, 
and § and k are the temperature gradient 
and the thermal conductivity respectively, we 
can deduce the relations, 
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Fig. 3. The records of the adverse temperature gradient obt 
The letters “on” and “off” indicate the duration of the heating 


heating currents were 0.500 amp. for a, 0.800 amp. for b and 1.323 amp. 


(95) 
(96) 


Tiix K 


and dnt=R Pz 


where #’ is a certain constant. 
The records of type (b) and (c) which reached 


“ ales £ 
a ett cect secant arte 


ST 
rigeiie Bibi 


? 
4 
: 


ained in the experiments of d=6cm. 
applied at the bottom. The 
for c, respectively. 


certain maximum values and decreased for the 
continuous heating indicated the transport of 
warm fluid to the upper part of the layer 


and they were classified as the records of 


convective heat transfer. The plotting of the 
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the bottom. 


maxima (8») of the temperature of the records 
(b) and (c) against the corresponding rates of 
heating in the similar way as in (a) give the 
segments like ## and vv in Fig. 4 respectively. 
The main difference between the records (b) 
and (c) was that records of type (c) showed a 
damped oscillatory trend while type (b) had 
no such oscillation. This difference of records 
suggested the difference in the mode of the 
convective motion, which also are seen in the 
relations between 8, and /,2 in the diagram 
as the segment vv is almost linear while 7? is 
a smooth curve. 

Considering the method of measurement of 
the adverse temperature gradient, the differ- 
ence were interpreted as follows. If the con- 
vection started simultaneously over the whole 
horizontal dimension of the layers, the peak 
value of the adverse temperature gradient 
should be recorded as soon as the warm fluid 
reaches the level of lower set of junctions. 
Then, a sharp decrease of the adverse temper- 
ature gradient should be recorded when this 
warm fluid reaches the upper set of junctions. 
Consequently, under such circumstances re- 
cords of a sharp peak like (c) should be ob- 
tained with the onset of convection. While if 
the convection start at a certain part of the 


The relation between the adverse temperature gradient and the rate of heating applied at 
The data are taken from the results of experiments of d=6cm. 


fluid it should be subsided by the surrounding 
cooler fluid and the resultant slow motion of 
the convection would produce a round-out 
peak like the record type (b). 

Literally understanding the condition of the 
criterion (54) and considering the small dimen- 
sion of the container and the non-uniformity 
of the magnetic field, the data of the adverse 
temperature gradient obtained from records 
type (b) were omitted for the determination 
of the critical temperature gradient. Thus, 
the critical adverse temperature gradients 
which entered the RAYLEIGH numbers are de- 
termined by the intersection of the extended 
linear trends of the segment like dd and vv, 
after applying linear least squares approxima- 
tions for each of the segments respectively. 

Another approval of the above procedure 
could be obtained from the results of theoreti- 
cal study of this problem under tilted exter- 
nal magnetic field (S. CHANDRASEKHAR 1954b), 
because under such circumstances smaller 
critical values of the RayLeIGH number have 
been obtained with the cells in the form of 
roll. Therefore, it is plausible that the seg- 
ment 7¢ is caused by the nonvertical charac- 
teristics of the magnetic field. 

The final results of the experiments are 
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Table 1. 
Depth Conductive (dd) Be=hy TZ +c, Convective (vv) Pp=hyI,2+ ey Q R 
oP ky C4 a) Kee Cy Ox ~ 108 1S 
3 2.878 0.004 0.002 11.945 0.353 0.006 1.22 “1.35, 
4 2.791 0.001 0.001 11.249 0.202 0.003 2.18 | 2.50 
5 2.590 0.000 0.001 8.698 0.120 0.001 3.40 | 3.81 
6 2.784 0.000 0.001 12.047 0.093 0.001 4.00 | 5.65 


a; denote the standard deviations from the assumed linear relations. 


per °C/cm), c; and o; (°C/cm). 


summarized in Table 1 in terms of the critical 


-RayLeigh number and magnetic number Q 


with the coefficients of linear relations obtained 
in the least square fittings. The small stand- 
ard deviation from the linear relations should 
be easily noticed in the table, which suggests 
the reality of linear relation especially in the 
regime wv as discussed in the next section. 
The average intensity of the magnetic field is 
assumed as 1,400 gauss for all cases, consider- 
ing the results of the mapping of the strength 
of the field. 

The boundary conditions for our experiments 


10° 
—————_— Both rigid 
<a One rigid ,One free 
——____—— Bon ifee 

io 


Log R 
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The units are ke; (watt/cm 


were those for the two rigid surfaces, because 
of the formation of a contaminated film at 
the top surface preventing any motion at 
such surface. The results of experiments 
quantitatively confirmed the theoretical pre- 
diction as in Fig. 5. The curves in Fig. 5 are 
drawn from the theoretical results obtained by 
CHANDRASEKHAR (1952a) and the solid circles 
indicate the results of experiments. 


$8. Heat Transfer by Cellular Convection 


The theoretically noticed stationary charac- 
teristics in the relation between the critical 


io 104 10° 


Log Q 


Fig. 5. The dependence of the critical value of the Rayleigh number on the magnetic number Q. 
The depths of the layer of mercury are indicated. 
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value of the RayLEIGH number and a® (A. 
PeLttew and R. V. SoutHweE.u 1940) have been 
observed in several experiments (R. J. Scamipt 
and O, A. Sanpers 1938, Y. Nakacawa and P. 
FRENzEN 1955 and others). By visible observa- 
tions, the author has confirmed that the cells 
appeared at first remain over a certain range 
of rates of heating beyond the point of inter- 
section, such as Borir in Fig. 4. For the range 
of rates of heating corresponding to the seg- 
ment vv in Fig. 4 we can safely assume that 
the convective heat transfer is mainly perform- 
ed by the similar size cellular convection which 
appeared first when the critical condition is 
exceeded. 


Then the difference of the peak values of 
temperature gradient (8») of the records (c) 
are able to be interpreted as the difference in 
magnitudes of the vertical velocity (w) and 
the temperature fluctuation (@) of the convec- 
tion. Relating the energy of the convection 
to the external energy supply, we could assume 

1 a2 2 

5 Ow Sole (97) 
where bar is used to denote the average over 
the cells. The amount of heat transfer (/’) 
by convection can be expressed as 


+1/2 eee (le 
Ke« \ C,pw6 dz « | C,pwdzZ , 


-1/2 -1/2 
(98) 
where the last relation is obtained 
from the proportional character- 
istics between w and 0 as shown 
in the equations of the theoretical 


io 
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studies. Consequently, we can 
combine (97) and (98) as 
K«In?. (99) | 


While from the results of experi- 
ments we can assume the relation 

TAR By (100) 
along the segment vv, then using 


the relation (99) and (100), we 


can write 
Box Ow w , 


(101) 


Now, if we assume that the in- 
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itial well defined oscillation in the record (c) is 
due to the vertical motion of the warm fluid, 
we could relate the period (r) of the oscillation 
to the vertical velocity in the following 
manner: 


(102) 


where d is the depth of layer. Then combining 
(102) and (101) or (102) and (100), we have 


TY Bp=const. or t/,=const. (103) 


The relation (103) is satisfied by experiments 
as illustrated in Fig. 6, where the data was 
taken from the results of similar experiments 
(technical reports of the author). In other 
words, the experimental proof of (103) supports 
the relation (100) and (101). Then, from (99) 
we can write 

=k By. (104) 
The relation (104) is the relation which has 
been especially used in the determination of 


the critical value of the RAYLEIGH number by 
some of the experimental investigators. 


$9. Summary and Remarks 


The development of the studies on the in- 
stability of the layer of fluid heated from 
below under various circumstances has provid- 


“NB, = Constant 


VBp (Arbitrary Units) 


Fig. 6. The empirical relation of rt and By The data are 


taken from the results of experiments of rotating layer of 
water (d=3cm, 2=10 rev/min). 
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ed valuable informations in relation to its ap- 
plications to the problems of astrophysics and 
geophysics. 

The applications to the meteorological pheno- 
mena, especially to the form of clouds have 
been studied by D. Brunt (1937, 1951) and 
others. The earlier difficulties on the direction 
of the convective circulation in meteorological 
applications, that is, the downwards motion at 
the center of cell in the air, have been removed 
by the study of Y. Naxacawa and P. FRENzEN 
(1955). Also, further application of the results 
of such studies to meteorological problems 
including the overstable convection has been 
proposed by the study of the above investiga- 
tors. Anexample of the oceanographic appli- 
cation is the formation of stripes on the sea 
surface which has been observed and explained 
by A. H. Woopcock and J. Wyman (1947) as 
the results of cellular convections over the sea. 

The experimental studies in cases of the 
electrically conducting fluids are especially 
significant because they have provided the 
justification on the study of the problems of 
magneto-hydrodynamics, and the only means 
of study on such subjects have been limited in 
theory. The applications, such as discussions 
on the problems of sunspot, magnetic stars 
and geomagnetism have been presented by S. 
CHANDRASEKHAR (1953d) and W. M. ELSAsSsER 
(1955, 1956). 
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The assumption of the linear temperature 
gradient in theories before the onset of con- 
vection has been justified by experimental 
investigations as that assumption was found 
to be satisfactory through the layer of fluids 
except in the very narrow region of the 
boundary of the layer (technical report of the 
author). 

A discussion on the possible shapes of the 
cell was presented by H. Stommen (1947) with 
the assumption of tessellation of plane by 
identical poligons. He showed that the shapes 
of equilateral-triangle, square and hexagon 
were the shapes which satisfied such assump- 
tion. In most of the experiments, the pre- 
dominance of the identical hexagon has been 
reported. The particular dimensional propor- 
tion of the cell in terms of a? has been 
experimentally confirmed in various circum- 
stances (see the results of R. J. Scumipt and 
O. A. Sanpers 1938, Y. Naxacawa and P. 
FRENZEN 1955). 

The effect of the motion of the fluid was 
studied theoretically by H. Jerrreys (1928) 
who showed that except the possible deforma- 
tion of cell in the direction of flow, the criterion 
would be the same. The effect of a shearing 
flow in a certain direction or the presence of 
horizontal temperature gradient has been studi- 
ed by some authors, and in some experiments 
cells were deformed to roller patterns in the 


Table 2. 
| “4° 
Exchange Stability Overstability 
Q=0 Re is a fixed value. | 
1 —() Experiments confirmed. | 
= | “ : R¢* « T?/3 (when TP) Experiments confirmed, 
Oe R.« T?/3 (when To) Experiments confirmed. g The validity depended upon P and ‘< 
> t A : = Re*«Q (when Q>) Only feassible for 
ie Rex @ (when Qo) Experiments confirmed. ‘ astrophysical conditions, 7<K. 
| — = — = a = — = 
Re«Q (when Q>o for fixed T) Re*«Q (when Qo for fixed 7) 
ave No experiment. _No experiment. 
Two different modes become possible The validity complicatedly depended 
T~0 depending upon 7 and Q. upon 7, P, Q, and w. 


Throughout the table the difference due to 


the boundary conditions are understood. Fe and 


R-* are used to denote the critical values of the Rayleigh number based on the exchange stability 


and overstability respectively, and other notations denote similar meani 


tical studies. 


ngs defined in the theore- 
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direction of flow (K. CHanpra 1938) or the 
direction of temperature gradient (Y. NakaGa- 
wa and P. Frenzen 1955). However, with the 
plausible explanation like the discussion on dé 
in section 7, we can easily interpret the essen- 
tial feature of the problem. 

Finally, the present state of the studies and 
the results of studies are summarized in 
Table 2. 
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Summary 


The kinetic theory of gases for the rotating system is investigated with an appropriate 


BOLTZMANN’s equation for such a system. 
The problems of diffusion, hydrodynamics and thermal 


ed in BOLTZMANN’s equation. 


Thus, the COLIORIS’ acceleration is includ- 


conduction are studied with the derivations of the coefficients of diffusion, viscosity, and 
thermal conduction from the microscopic view point. 


Introduction 


The macroscopic properties of gases and the 
governing equations of such properties have 
been interpreted from the microscopic view 
point. By virtue of the kinetic theory of 
gases (for examples Sir J. Jeans 1940) and 
liquids (M. Born and H. S. Green 1949), these 
macroscopic properties are shown as the 
“averaged ’’ properties of molecules, also the 
equations are explained as the “‘averaged” 
movements of such properties. 

In many fields of astrophysics and geophy- 
sics, the rotating coordinate system has _ be- 
come essential for the frame of reference for 
_the study of the motion of such “‘averaged’”’ 
molecular properties. Some investigators of 
the kinetic theory of gases have paid atten- 
tion to the problems of the rotating gases 
(see S. CHapman and T. G. Cowrine 1939). 
Their treatments, however, were only referred 
to the fixed system and their interests were 
confined to the additional external force po- 
tential due to the centrifugal acceleration 
which came into the velocity distribution 
function for the equilibrium state. 

In this paper, the detailed study of the kine- 
tic theory of gases for the rotating system is 
undertaken and the appropriate representations 
of the ‘‘averaged’’ properties and the govern- 
ing equations of such properties are investi- 
gated for such a system. 


The BOLTZMANN’s Equation for the 
Rotating System 


Sale 


The BontzmMann’s equation for the rotating 
system can be derived in the following manner. 
Consider the number of molecules whose 
velocity is within the range c, e+de and the 
position is within the range r, r+dr at timez 
in the velocity-position space and let the num- 
ber be 

fice, r, Ddedr. 
After time d¢, this number of molecules would 
form the set 
fietedt, r+rdt, t+dt)dedr 

due to the motion of molecules, where the 
volume occupied by the set would remain the 
same as dedr by the Litouvinun’s theorem. 
However, during dt, some molecules would 
be deflected out of this set by collision with 
other molecules, while some other molecules 
would enter this set as results of collisions. 
This change should be proportional to de dr dz. 
Hence if we denote influences of collision by 
Of dedrdt, then, as the result of the net 


Ot 
balance of the number of the molecules in the 


specified set during the time dt, we have the 
following relation 


Yel de dr dt = { fle+ edt, r+rdt, t+dt) 


0 
—fiec, r, t)}de dr dt 


Offa aOfae. Op 
wt Sl ded 
case apt Fg) }dedr dt 
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where the following relations are used to 
obtain the final equation, 
T=Ch (1.2) 
c= F+2(c/\Q) , (1.3) 
and F' denotes the force depending on the 
position, and Q=(@,, @y,, @:) is the rotation 
vector. 

It is immediately apparent that the term 
corresponding to the Cortouis’ force becomes 
involved in the BoL_tzMann equation. The 
change in the BontzmManN' equation can 
be better understood when we consider the 
trajectory of molecule, as it is no longer a 
straight line, but becomes a diverging spiral 
for the rotating system. 


§2. The Motion of Molecules and the Dyna- 
mics of Collision 


Consider a typical rotating system such that 
the coordinates rotate around z-axis with a 
constant angular velocity w**, then the equa- 
tions of motion of a molecule in such a system 
are, 


mx=2mwoy+ moe , (221) 
my = —2mox% +My , (2.2) 
mz=0, (2:3) 


where coordinates are taken in the rotating 
frame and the external force is omitted in 
order to see the essential features of the 
results of rotation. The solutions of the equa- 
tions (2.1), (2.2), and (2.3) are easily obtained 
and with initial conditions at 7=0, 


c=, Y=Yo, 2 =o, il 


1h See (2.4) 
=2%0, Y=Yo, Z=2p, J 
the solutions are 
X=x COS ot + yy) Sin wt 
+4{(xp— WY) COS Wt + (Yo— 2x) Sin wt}t , 
(2.5) 


* Throughout the paper the notations are used 


after S. CHAPMAN and T. G. COWLING 1939. 

** Through the paper such a coordinate system 
is mostly used. However, the results obtained can 
be easily generalized. 
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Y=Yo COS wt—29 Sin wt 


+{(yg—w2q) cos wt —(%p— @Yo) SiN ot}t, 
(2.6) 

Z=2y thot. (2.7) 
Thus, for diverging spiral, with constant veloci- 
ty 2) along the z-axis, these give a planer pro- 
jection of a diverging spiral on the zy-plane. 
(see Fig. 1) The solutions (2.5) and (2.6) 
suggest that the phenomena observed in the 
rotating coordinates would be somewhat simi- 
lar to the manner in which ionized gases _be- 
have under an external magnetic field.*** In . 
eqations (2.1) and (2.2), if we neglect the last 
terms in the equations, the solution would 
give a closed spiral as it is obtained in the 
ionized gases. ; 

The change in the trajectory brings up the 
problem of the possible change in the dyna- 
mics of the collision. However, in the 
dynamics of the collision, we can assume 
that the time duration of collision as infinite- 
simally small. Then the equations (2.5), (2.6) 


Z 


x 
Fig. 1. The trajectory of atypical molecule in 
rotating coordinates. P denotes the initial 
position at t=0. 


eK see ibid. Ss CHAPMAN and TG COWLING 1939, 


chapter 13. 


. fixed system. 
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and (2.7) would be approximated in the fol- 
lowing manner: 


Y=Yot yott+OW) , (2.9) 
2=2)t2ot . (2.10) 


From these equations, a linear motion of a 
molecule during the collision would be able 
to apply for the consideration of the dynamics 
of the collision. Consequently, the dynamics 
of the collision be the same as that for the 
Hence, the number of mole- 
cules, linear momentum, and the total energy 
would be conserved during collision. With 
the conventional notations of the function ¢™, 
we could write these collision invariants as 
below: 


the number of molecules; 


pO=1 (2.11) 
the linear momentum; 

g@=mC, (2.12) 
the total energy; 

Oey ; mC, (2.13) 


where C, C denote the vector and the scalar 
expressions of the peculiar velocity of the 


molecule defined by 
C=c—o, (2.14) 
_ where 
ney=\ fe de, (2.15) 
and ec denotes the ‘‘averaged’’ peculiar 


velocity, that is, velocity of the mass motion 
of gas, and a bar means that the quantity is 
‘*averaged’’ over the whole velocity space. 
Rewriting the BoLtTzMAaNN equation (1.1) in 
terms of the peculiar velocity, we have 


1OyF of ca], OF. 
Or Gr +| F+2(6/\@) - sea ac 
of of Of (916 
+2(C/\2)- 56 — 56 ‘gy “Fer hae, 
where 
D 0 0 
DE hor 


is called the ‘‘ mobile-operator ”’ 
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A similar equation can be obtained for a 


function of molecular properties ¢, after 
averaging over the whole velocity space, 
Dno (Db Op 
ene +2. = 
Dt ap Ot mec neon Or 
Dey | 06° 
F+2 Q)— ~~" |. 
+ +2(¢9/\ Q) BoA ac 
O08 OP, LUGAR) ON 
OC Or eal 
=| Dfde, (2.17) 
where 
‘ goed 
G = 
| 2 fae=\* OL de 


and the bars denote the averaged values over 
the world velocity space. 

Substitution of the encounter invariants for 
the function ¢ gives special forms of the equa- 
tion (2.17), which are the fundamental equa- 
tions in physical problems. 

The substitution 


p=$O=1 
results in the equation of continuity, 

Dn 0 

N—-Cy= Pe 

Dt cs Lie oar an) 
or 

D 

a +05 -c0=0, (2.19) 
where 

p=mn 


denotes the density of the gas. 
The substitution 


o=$=mC 
results in the hydrodynamical equation for 
the rotating system, 


Deo_ eee 
2(eo/\2) 5 ar 


iB -P+F, (2.20) 
Dt 


where 

P=pCC 
is the pressure tensor. This equation is 
formally identical with the result obtained by 


the ordinary transformation of coordinates, but 
as is shown later, the quantity represented 
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by the pressure tensor is not immediately 
identical with the result derived by the co- 
ordinate transformation. 

The substitution 


o=~@O=E 
results in the San of thermal conduction 
DT Zot 0 
— = —- ; ae 
Dt 3nk bas ee a} ( ) 
where 
k=” Ga? T 
2 9 
and q=nCE 


denotes the energy flux vector. Again the 
resultant equation bears similarity to the ex- 
pression derived for the fixed system. How- 
ever, the pressure tensor P and the energy 
flux vector gq would become different from 
those for the fixed system as are shown in 
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later sections. 
Equation of angular momentum conservation 
can also be obtained, when we _ substitute 


mr /\c for ¢. 


§3. The Pressure Tensor and the Hydro- 
dynamical Equation 


The BoutzMANN equation (2.16) is solved by 
the method of successive approximations deve- 
loped by Ensxoe and others, and the coefficients 
of thermal conduction, viscosity, and diffusion 
are obtained in this process of approximation 
in the similar manner for the ionized gases 
studied by T. G. Cownine (1932, 1933) to show 
the effect of rotation. 

Suppose the temperature is uniform through- 
out the gas, the components of the pressure 
tensor can be obtained in the following forms 
after a little laborious calculation 


ie} 
P= a mt = | e205 + Gert) 8 (4e02 Vee *s ree 3 Gar |, (3.1) 
+ EZ w? v2 
9 
&) Que 2 Ih 2 
Pyy a 16 (4 ' Ez a= D) Gx +8 10) 6 4652 \c +b ry cs 3 20 | A (3.2) 
+ w? v2 
9 
Du ena (3.3) 
° ° ° 
Pay =pygD = — — ee gE ue 1 (Gx —Eyy) é oye | P (3.4) 
+ 9 (4w?)¢? 2 3 
Po pO ea mens gE _2 Qos. | : (3.5) 
5 (4 w?) tc? 3 
° 
Dy mea Pet ae bs — je. + : (2a) ee. p (3.6) 


ap ; (4?) t? 


where 
_ 3 
Ty a 
° ° : 
and @xz, @yy-+-etc. denote the components of 
fo) 
tensor (0/Or)e,, and the rotation vector 


2=(0,0,) is taken in z direction. The super- 
script (1) denotes the order of the approxi- 
mation, namely in this case the expressions 
were obtained up to the first approximation 
of the pressure tensor. y is the viscosity 
being defined to the first approximation in 
ordinary theory, 


In. the above expressions ‘of the pressure 
tensor, all the components except p.2™, indi- 
cate the influence of the rotation of the co- 
ordinates. However, under the usual circums- 
tances the factor r is . 


310" 
2x 10° 


and negligible in these expressions. With 
this order in mind, the hydrodynamical equa- 
tions for the rotating coordinates which is 
rotating around the z-axis with a constant 
angular velocity » should be written in the 


~ O(10-?) eae 
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following form 


Du ibe: a) 

DE ee = ce an 4 Ds. 8+ vy ety to? wt+F,, 
(3.8) 

Dy 10 

pt teem == baat : * O+vpn+toyt+F, , 
(3.9) 

Dw Wo. 4) 

a am ate 3 Oe -O+yvp?wotF:, (3.10) 

where 

p=nkT 


is static pressure, 
\3 ie 


p 


is the kinematic viscosity, and the notations 
9 and p? are 


Ou Be Ovo 4 OWo 


@= 
Ca Git Og? 

ke 2 — ids 
G2? Oe Oe" 


Hence, the hydrodynamical equations deriv- 
ed by the transformation of coordinates are 
practically secured for its use in the investiga- 
tion of the problems under ordinary circums- 
tances by virtue of the foregoing discussions. 


$4. Thermal Conduction 


The thermal conduction for the rotating 
system is considered in the following manner. 
Suppose the gas is at rest relative to the ro- 

_tating frame, that is, the gas is in a state of 
rigid-body rotation. 

If the temperature gradient 07/0r is the 
direction of the rotation vector 9, it can 
easily be’ seen in analogy to the ionized gas 

that the heat conduction parallel to the axis 
of the rotation is not affected by the rotation 
of the coordinates. 

In the other special case where the tempera- 
ture gradient is perpendicular to the vector 
of rotation, we have for the expression of 
the enegry flux vector g™ 


Os 
if ar —2r (25, ar =a ip 
1+ 4a?r? 


where superscript again denotes the expression 
held to its first approximation and 


qg?= = 


Own OL 0 f 
or =( Ox’ Oy ’ 0) 2= (0, 0, w) 
and 
at . 
2p 


A is the first approximation to the ordinary 
thermal conductivity in the fixed system. In 
this case, the thermal conduction along the 
gradient of temperature is reduced to 
1/(1+4@?c?) times, while the other conduction, 
the ‘‘transverse conduction’’, which is per- 
pendicular to both of the vectors 0T/0r and 
2, becomes in the magnitude 2wr times the 
former conduction as the results of the rotation 
of the coordinates. 
The order of the quantity rt is, again, 


34 34 ota 
=5g> ae Od0-*) (4.2) 
since 
i= Coie 


Consequently by the similar reason as in 
the previous section, the final form of the 
equation of thermal conduction under the pre- 
sent assumption would become 


DT 
— tf 
Di KV (4.3) 
where 
year 
3kn 
is the thermometric conductivity. 

However, if rt is large, the thermal conduc- 
tion would mostly be directed along z-axis. 
Back to the equation (2.21), the term (0/0r)-q 
are 


04x A Geaig, VanwAas oar 
= — ee 
Ox 1+4?r? Ox? 1440? Oxdy 
O”gy “thee A or 2cwd eT ‘ (4.5) 
Oy 1+4w?r? Oy? 1+402r? Oxdy 
odie ee (4.6) 
dz 422 


By addition of these three equations, the 
second terms in (4.4) and (4.5) would cancel 
each other. Therefore, the resultant gives 
such a conduction that heat would be conduct- 
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ed mostly along the z-axis, because the hori- 
zontal conduction becomes smaller by the 
factor 1/(1+4w?r?), while no transverse con- 
duction appears. 


$5. Diffusion 


Similar results as in the foregoing sections 
can be obtained in the problem of diffusion by 
considering a gas-mixture. If we take the 
case of binary-mixture, we would have the 
following equations for conservation of mole- 
cules, hydrodynamics, and thermal conduction, 


4 0 
OP ae neat nC + tn) =0yy, 0) 
Do 0 
-Co=0 5.2 
Dt CO a Co ( ) 


(The equation of conservation of mass), 


Dey PEA ano, 3 
P DE = 01F'+ 0.F 2+ 20(e/\ 8) Or Peds oes) 
3 in = S kn T 2 -(miG+ mC) 


+ pF: C,+ poke: C, 
6) 0 

ar Co yr , 
where notations are used in the usual way 
and the subscripts 1 and 2 refer to the lst 
and 2nd kinds of gases, respectively. From 
them, and with the assumption of solid rota- 
tion and uniform temperature of gases, the 
velocity of diffusion is obtained as 


_ pet dy,—2t(Q/\di) 
0102 1+40?? : 
This result indicates that the direct diffu- 

sion is reduced in the ratio 1/(1+4w?r?), while 
the transverse diffusion is 2r times the direct 
diffusion. Therefore, a result similar to that 
discussed in the last part of section 4 also 
holds in diffusion provided that 4w?r? is large 
to be comparable with unity. Under the last 
circumstance diffusion would be mostly direct- 
ed along z-axis. 


=p: (5.4) 


G-G= (5.5) 


$6. Summary 


The proceeding study has showed that the 
various dynamical equations in rotating system 
can be easily derived as the ‘‘averaged’”’ pro- 
perties of molecules of rotating gases, also it 
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confirmed the negligible influence of the rota- 
tion in microscopic characteristics of the 
gases, such as the coefficients of viscosity, 
thermal conduction and diffusion. The latter 
indicates, under ordinary circumstances, that 
the MAxweE.LuiAN velocity-distribution function 
should hold in the rotating gas, as well as 
the equipartition law of energy flow. Conse- 
quently, it becomes clear that the cause of 
the so-called two-dimensionality of the motion 
in the rotating gas and liquid must be a 
phenomenological one. 

This proof can be easily obtained as follows: 
If the approximation of the hydrodynamical * 
equation, 


o\a=—-12 p (6.1) 
00r 

holds with the condition of the incompres- 

sibility. Then by taking curl of this equation, 

we can obtain the following relation 


0 
(2-5, )au=0 . 


The equation (6.2) denotes that the motion 
only takes place in planes perpendicular to 
the axis of rotation, in other words, the 
motion has the characteristics of the two- 
dimensionality. 

In the large scale circulation of the atmos- 
phere and ocean, these conditions were found 
in close approximation, also in some experi- 
mental studies of the rotating liquids by G. I. 
Taytor (1921, 1923). However, the im- 
portance of the study of three-dimensional 
motion in the rotating system should be easi- 
ly understood from the results of the present 
study because the influence of the rotation is 
negligible in microscopic characteristics. 


(6.2) 
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Wave Groups Generated by a Very Small Explosion. 


By 


Kyozi Tazime 


Department of Geophysics, Faculty of Science, 
Hokkaido University, Sapporo. 


Summary 


Generation mechanisms of wave groups by a small explosion are studied from a 


wave theoretical point of view. 


We see all waves predominant in the seismic records 


obtained are some kind of surface waves, especially of dispersive RAYLEIGH and SEZAWA 
waves. These wave groups are shown to satisfy the quarter wave length law which has 


been established theoretically in the present paper. 


Owing to shielding effects of a super- 


ficial layer, not very many kinds of surface wave groups can grow. This makes our seismo- 


grams rather simple. 


In seismic prospectings, the surface waves studied in the present paper may be useful, 


at least, in some cases. 


Introduction 


§1. 

When geophones are set up on a straight 
line at equal intervals in field experiments, 
we observe that all peaks and troughs of every 
trace in seismic records correspond very well 
to each other. These peaks and troughs will 
result in the respective time-distance lines 
which have one constant slope for one group, 
but another constant slope for another group. 
In general there are only a few groups having 
equally sloped lines. 

Within a wave group, we can observe 
several peaks and troughs at almost equal 
time intervals. As we have a shock type 
pulse at the origin, we are apt to attribute 
the generations of these oscillations in a 
group to dispersions of seismic waves. In 
this case, however, we soon meet with the 
difficulty of explaining the constancy of period 
in the wave group. 

As we have got more or less similar seis- 
mic records in our experiments to be described, 
we shall confine ourselves to the study of 
records obtained in 1953 at the ground of 
Geological Survey of Japan. 


§2. Seismic records 


At first a cap was blasted at the depth of 
0.75m. Geophones were set up on the ground 


at every 0.25m (or 0.50 m) on a straight line. 
The epicentral distance spreads from 0.25m 
fo< 25:25 eb. Le CE) ES. Se Gz Oeand 
Sassa-SANEl’s (U) seismic recorders were used. 

Seismograms thus obtained are shown in 
Fig. 1 and Fig. 2. Since various kinds of 
recorders were used at the same time in this 
experiment, the forms of these seismograms 
appear to be different from each other. The 
wave forms are, however, out of our concern 
at present. 

Fig. 3 shows the time-distance curves for all 
troughs contained in the above seismograms. 
Those for the peaks are omitted to avoid a 
confused figure. This chart is, the author 
believes, useful for the classification of wave 
groups. The slopes of the lines connecting 
each trough indicate the phase velocities of 
the respective waves. Thus we have four 
wave groups having different phase velocities 
as in Table 1. 

We have stated about the vertical component 
of ground motions alone, but the horizontal 
components were also observed. These are 
shown in Fig. 4 in which the transversal 
component is seen very small in all wave 
groups. On the other hand, the longitudinal 
component is predominant in the 3rd and the 
4th groups, though it is also poor in the lst 


ee oo 
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Fig. 1. A cap was blasted at the depth of 0.75 m. 


The recorders G and O were used. Geophones 


of vertical components were set up from 0.25m 
LOM Oamiae 
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Fig. 2. Vertical U geophones were set up from 
1 25em tor2on2o me 


and the 2nd groups. 


The Ist group represents initial motions and 
contains a few oscillations even at distant 
observing points. The other groups growing 
gradually with the epicentral distance, on the 
other hand, show several oscillations. The 
wave form of the lst group is very simple 
and this initial motion is interpretted to be 
due to the direct or refracted P wave. From 
now on, we shall concern with the wave 
groups except the first one. 

The 4th group is conspicuous everywhere 
including the epicenter and seems, at first 
sight, to be due to a kind of RAYLEIGH type 
surface waves, though its oscillatory character 
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Fig. 3. The time-distance curves of all troughs in 
Fig. 1 and Fig. 2. 


Table 1. Classification of the wave groups. 
group phase velocity period 
I | 1.3~18 x 10? m/s | 
II | 5.0~6.5 2.8 x 10-2 sec 
Se eee 2 6 2.5 
IV 0.6~0.7 | Be) 


| 


has not yet been thoroughly understood. 

It is more difficult to understand the genera- 
tion mechanisms of the remaining two wave 
groups, the 2nd and the 3rd. The phase 
velocity of 2.5x10?m/s, for instance in Fig. 
3, is observed within 6m from the epicenter 
as well as beyond 11m. But it cannot be 
decided from Fig, 3 alone whether the both 
tremors belong to the same wave system or 


15) 


Fig. 4. The upper two; 7, L 15.75~17.25m, 
and V 15.75 ~ 18.75 m. 
The lower two; 7, 017.25 ~18.75m, and V 17.25 


~ 20.25 m. 


not. The same phase velocity in the two 
parts may be an accidental coincidence. 

It must be remarked that the word “‘ wave 
group ’’ in this paper has no complete physical 
meaning, but merely the group 
having the proper phase velocity. The full 
discussions on this point will be given later. 


indicates 
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§3. Rise and Fall of the Individual Peak 
or Trough within One Wave Group with 
Regard to the Epicentral Distance 


Suzuxi and Sima (1954) found an apparent 
rule as to the form of each phase recorded on 
just the same seismogram as in Fig. 2. 
Modifying the Suzux1 and Stma’s method, we 
have tried to find the feature of rise and fall 
of the individual peak or trough within one 
wave group with regard to the epicentral dis- 
tance. In order to do this, we read ampli- 
tudes of all peaks and troughs in one trace of 
the seismic record. If fortunately wave groups 
are separated form one another as in Fig. 2, 
we take the sum of all amplitudes of peaks 
and troughs in a group as an amplitude 
standard of the group. The ratio of each 
amplitude to the sum is a standardised am- 
plitude of the respective peak or trough with- 
in the wave group. Other trace of the 
seismogram is treated in the same manner 


S aoe j 
WD Sta at a ee 9d") tia 3u + 2,0 
| @u F on sesettag 00 
| ed : eS ee eae 3d ee 
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Fig. 5. The amplitudes of peaks and troughs 
within wave groups. 
The abscissa indicates the epicentral distance. 
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and so on. In doing so, the difference in 
amplification characteristics of recorders may 
be eliminated at least in the first approxima- 
tion. Thus we can follow the rise and fall 
of every phase with respect to the epicentral 
distance. The other wave groups can be 
analysed similarly. Fig. 5 illustrates the fea- 
tures of rise and fall for the 2nd, the 3rd and 
the 4th wave groups. The notation wz or d 
in the figure corresponds to up or down on 
the seismograms and the number corresponds 
to the order of the peak or the trough from 
the initial motion. : 
We see in Fig. 5 that one early appearing 
wave grows gradually at first up to a certain 
epicentral distance, beyond which it decreases 
and the next wave comes to be predominant. 
This feature is seen everywhere from the 
2nd group to the 4th group. In order to see 
this more clearly, we combined the two 
figures, Fig. 3 and Fig. 5, and obtained Fig. 6 
in which the above mentioned distances are 
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Oo 
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Fig. 6. The two figures, Fig. 3 and Fig. 5, are 
combined. 
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marked by small black circles on the time 
distance lines. 

We see that the origin of the wave having 
2.5 10?m/s phase velocity can be traced back 
to very near the epicenter. The broken lines 
in Fig. 6 indicate the parts where we cannot 
trace the wave, because of overlapping of 
wave groups. The meaning of the chain lines 
will be explained in the next section. 


$4. The Group Velocities Directly Derived 
from the Seismograms 

Looking at the form of a wave group on 

' our seismogram, for instance in Fig. 2, we 

see it bounded by a spindle-shaped envelope 


aseim Ho. #7. 
; Te 
Dt 
a ee 
Fig. 7. Each wave group is enveloped in a 


spindle-shape. 


Tracing the loops of these spindle-shapes, 
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Fig. 8. The amplitude of peak or trough with- 
in a wave group changes with respect to 
time as well as epicentral distance. 


Thus we see the maximum amplitude of a 
wave group is propagated with a certain defi- 
nite velocity. Thus the slopes of the chain 
lines show the minimum group velocities of 
the wave groups. The group velocities Us, 
thus decided, are tabulated in Table 2. 

For reference the frequency-responses of 
our recorders are shown in Fig. 9. 


we get the chain lines drawn in Fig. 6, Table 2. Observed group velocities. 
though the starting points of these lines are group | group velocity (Uo) 
somewhat obscure. At any rate, it must 2 eee = 
be noted that the chain lines pass through i 1.310? m/s 
the black circles already described. This I 0.92 
situation is shown schematically in Fig. 8 IV 0.52 
lu 
Ss |.0 x a 4 
= = 
: a a3 
s = 45) = = 
5 Ww ies u 2 
eG 2 > 
‘Ss 1000 BE E 
ob a ie <q 
oO a | 
& 500 we [ae eee et 
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Fig. 9. Frequency-responses of our recorders. 


§5. Analysis of the Wave Forms 


Notwithstanding hardly noticeable change 
in periods within the wave groups, dispersion 
phenomena will be studied here. If we assume 
that the spindle-shaped wave form as in Bice: 


is composed of two sinusoidal waves having 
equal amplitudes, we have 
y=cos (2z#t/t)-+c0s (2n¢/t") 
—2 cos (2nt/r) cos {nt/(ce’/|t—7’|)} , 
where ¢ means time and r and c’ are periods 
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which are nearly equal. Therefore the period 
of an envelope of the waves is 


T=2re'/|t—t'l 

Now, for instance, for the 2nd wave group 
in our seismograms, we have 

Geno eA SEC 
and 

(= UPN Se. 
hence 

C =Tit e ot) =2.0 Of a4 x LO Sec... 

Thus we see the wave group has a very 
narrow period range outside of which the 
amplitudes of individual waves will fall off 
sharply. In practice, it will not be possible to 
distinguish the change of periods between 
25x10 sec. and 2.8x10-7sec on an actual 
seismic record, Thus we must give up to 
expect dispersion effects as in usual seismo- 
grams by natural earthquakes, since we have 
no changing period within a wave group. 

In order to make this point clearer, FourIER 
analyses of wave forms were carried out. 
The result of the analyses for the 2nd wave 
group is shown in Fig. 10. 


o 27t=57xl0%Sec 
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Fig. 10. Result of the FOURIER analysis. ©, A 
and x are calculated from the same seismogram, 
but the time units used are different from one 
another. 


In this figure 2x means the time interval 
adopted for the analyses. We have more- 
over tried to find phase velocities of individual 
FourIeER components. But this attempt is not 
successful, since enormous errors enter in the 
culculations of phase velocities outside the 
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above mentioned narrow period range. It 
seems, however, that the predominant phases 
in our wave groups have the phase velocities 
near the inflection points on the dispersion 
curves. Estimating the group velocities with- 
in the narrow range of period from the cal- 
culated phase velocities, we obtain respective- 
ly almost the same values as those in Table 
2 which were directly derived from the seismo- 
grams. This discovery makes us_ believe 
that the tabulated group velocities are the 
minimum ones, as are predicted by theoreti- 
cal studies of this kind of problem. 


$6. The Subsurface Condition of Our Ex- 
perimental Ground 

By now we have investigated only the wave 
forms appearing in the seismic records and 
have not yet discussed generation mechanisms 
of the wave groups. Now the subsurface 
condition of our experimental field will be 
described. We shall start with the discussion 
of the results obtained by refraction shootings, 
though there are several boring data as 
shown in Fig. 11. 
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Fig. 11. The boring data at the experimental 
field. 
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_ A reverse shot was done with a cap keep- 
ing the other conditions, such as shot depth, 
positions of geophones and so on, the same 
as before. The time-distance curves for the 


reverse shot, shown in Fig. 12, were found 


quite similar to those by the normal shot in 
Pig. 3. 
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Fig. 12. The time-distance curves for a reverse 
shot. 


Therefore we may regard the subsurface 
condition to be uniform horizontally. Then 
the underground structure may easily be de- 


119 


duced from the travel time of the initial 
motion. The result obtained is given in Fig. 13 
(a), where the numerals indicate the velocities 
of P waves in the respective layers, the unit 
being 102 m/s. Another result may be obtained 
if we assume a linearly increasing velocity in 
the superficial layer. This is shown in Fig. 
13 (b). The other result (c) was also obtained 
by Hatakeyama (1954), one of the member of 


our experimental group, for a deeper shot at 
the same ground. 
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Fig. 13. The subsurface seismic models. - 


Now we encounter with a difficulty in choos- 
ing the model to be adopted. Hatrakryama’s 
model gives the time-distance curve shown by 
the broken lines in Fig. 14, assuming a surface 
shot. On the other hand, the full and chain 
lines correspond respectively to the models (a) 
and (b) in Fig. 13, small black circles indicat- 
ing the observed values for initial motions by 
a near surface shot. In spite of rather large 
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Fig. 14. Time-distance curves calculated for 


various subsurface models. 


differences among the subsurface models, 
there are little differences among the time- 
distance curves obtained. As is easily seen, 
this is unavoidable as long as the initial time- 
distance curve alone is used for the subsurface 
analysis. At any rate, we cannot choose one 
among the three models given in Fig. 13. 
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As the explosion is very small, S-phase ap- 
pears scarecely in our records. Thus we 
have no data about S-phase. Fortunately, 
however, Muraucut (1955) has observed, by 
his special technique, the velocity of S wave 
in the superficial mud layer, and got Poisson’s 
ratio (c) from 0.21 to 0.24. As are seen in 
Fig. 15, the values of v,/vs and vpr/vs are 
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Fig. 15. Relation between vp, vs, vr and 


POISSON’s ratio for a perfectly elastic 
medium. 


nearly constant within the range of Porsson’s 
ratio obtained by Muraucui. From now on, 
we shall assume Porsson’s ratio to be 0.25 in 
our calculations. 

Thus the velocities of S and RayLEIGH waves 
can be estimated by using that of P wave 
given in Fig. 13. The results are tabulated 
in Table 3. 


Table 3. Velocities of P, S and RAYLEIGH 
waves for various models. 


mee an. (c) 
layer |P |S R|P|S|R|P|S|R 
Ist | 1.3 0.8) 0.7, 1.3 0.8 | 1.8 0.9 0.8 
2nd | 2.5) 1.2| 1.1) 8.0 4.7 8.8 5.2) 4.7 
3rd_— 18 u O58 11 | OSS iL) eb 

Table 4. Data on later phases. 
wave] component. [period Phase, [wave] group. 
iT | vertionl | 2.8 6.0 )i7mh ae 
I Vencnuen | 2 ole 2.2 | 5m 0.9 
v |verticll =| 35] 0.7 | ami 0.6 


| longitudinal 


time unit 10-2sec; velocity unit 102 m/s 
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Some data on later phases are collected in 
Table 4. All of the later phases have very 
small phase and group velocities, exhibiting 
remarkable dispersions. 


§7. The Quarter Wave Length Law Applied 
to the Superficial Layer 


We shall now compare records obtained when 
caps were blasted at various depths. Fig. 16 
was obtained by a shot at depth 1.50 m, the 
other conditions being kept the same as in 
Fig. 4 (d). We see the 3rd and the 4th 
groups become very poor while the 2nd group ‘ 
maintains its predominance. These observed 
facts are summarised in Table 5. Considering 


2 i Lie és 


Fig. 16. A cap was blasted at the depth of 1.50 m. 


Table 5. Relation between the predominances of 
the wave groups and the shot depth. 


shot depth (m) 


wave 
ia 0.08 OS: .- aes 
| 
Il F FO Ueag 
Ill F F P 
IV | F = oe P 


fair SPs poor 


these results in connection with the sub- 
surface models (Fig. 13), we find that some 
discontinuity plane exists between the depths 
of 0.75m and 1.50m. Thus we must give up 
the model (b) in Fig. 13. We know also that 
the 3rd and the 4th wave groups are related 
almost entirely to the superficial layer, and 
on the other hand the 2nd group is related to 
the lower layers. As transversal components 
are very poor for every wave group, all 
wave groups we have seen are considered to 
be dispersive RAYLEIGH or SEzawa waves, 
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usually called M-waves. 

In the case of a single layer on a semi- 
infinite medium, fortunately, several disper- 
sion curves of M-waves have already been 
calculated by Sezawa (1927), Kanai (1951) and 
others. As Toutsroy and Evucrene (1953) ex- 
plained, there are two groups of waves called 
M, and M,, and an additional suffix is used 
to denote the order in the respective groups. 
For instance, My, or M.; means respectively 
Mi or M,-wave of the first order. Now 
several dispersion curves have been reproduc- 
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waves in our seismic records seemed to have 
the minimum group velocities for the respec- 
tive wave groups. In Fig. 18, are shown the 
period, phase and group velocities of the Airy 
phases as functions of the rigidity ratio /2//41. 
We see that each curve in the uppermost 
figure has its asymptote as /2//4,;->©, that is, 
the period of each Atry phase becomes con- 
stant. The prefixes of M mean respectively 
the first and the second minimum group 
velocities, and Z,; corresponds to the first 
order Love waves. 

The asymptotic values of the Arry phase 
periods, ZT) say, are 
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ed in Fig. 17 where, unlike usual ones, 
log-log scales are used and abscissae are con- 
verted from wave length to period having a 
time unit H/vs. H means the thickness of 
the superficial layer and vs; means the veloci- 
ty of S-wave in this stratum. Log-log scales 
will make the comparison of calculated dis- 
persion curves with those observed easier. 
The conversion from wave length to period 
will be useful since we observe on a seismic 
record the period and not the wave length. 
Here we must recall that the predominant 


(b) Ma 
Theoretical dispersion curves of dispersive RAYLEIGH waves (Mi, M2:), 
by SEZAWA and KANAI. 


To/(H/ vs) =(L/H)(c/vs1) = 2.4 for Mi and Mp1, 
=0.8 for 1Mo1, 
=a) for UB 


where Z means the wave length concerned. 
We have previously assumed in Fig. 18 that 
Porsson’s ratio is 0.25, namely 
Un= Ws F3 

Therefore the above relations may be rewrit- 
ten as 

H=(1/4)vp1T, for Mu and .Mnp,, 

H=(3/4)vmTo for M21, 


We7? 
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Fig.18. The period, phase and group velocities 
of the AIRY phases against the rigidity 
ratio p/p. 


He WAVE T 5 for Vey, 4 
These relations should be understood as quarter 
wave length laws for doubly stratified layers. 


§8. Discussions of the Later Phases 


Putting the observed values vp, and T, in 
Fig. 13 and Table 4 into the above relations, 
we have 


0.8m for the 3rd group , 


i= 
; 1.1m for the 4th group. 


These are well coincident with the thickness 
of the superficial layer deduced by the re- 
fraction method, the latter being 0.9m to 
1.5m as in Fig. 13 (a) and (c). Comparing 
further the lower two figures in Fig. 18 with 
the observed data in Table 4, we see that the 
observed phase and group velocities for the 
3rd and the 4th groups are also well corres- 
ponding to the theoretical values in Fig. 18, 
if we take the rigidity ratio 4 to 5. Thus 
we see that the 4th group is the Airy phase 
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of My-wave and the 3rd group the second 


minimum of M,,-wave in the superficial layer 


in Fig. 13. 


In the above discussion, we have assumed — 


the rigidity ratio of about 4. This rigidity 
ratio seems adequate, since we get Up:/Vpi of 
about 2 by our assumption on Porsson’s ratio, 
no density jump being assumed. Therefore 


the subsurface model (a) alone in Fig. 13 will 


survive from now on. 
In the above, we have been able to explain 
the 3rd and the 4th wave groups by the two 


layer model, though actually there must be* 
But the observed — 


three or more layers. 
phases of the 3rd and the 4th groups appear 
to be scarcely influenced by the existence of 
lower layers, since the wave lengths of the 
phases are as short as 4 to 5 meters, being 
the same order as the depth to the third 
layer. 

The 2nd wave group has the more persistent 
period and the more regular wave form than 
the 3rd and the 4th groups. Thus the 2nd 
group seems again to be the Arry phase of a 
dispersive RAYLEIGH or SEzAWA wave. As the 
wave length of the 2nd wave group are longer 
than those for other wave groups in Table 4, 
we may expect the 2nd group will penetrate 
to the second layer, as shown by the full 
line in Fig. 19, whereas as the broken line 
shows, the 3rd group, for instance, will not. 
We must also keep in mind the fact pointed 
out in Table 4 that the longitudinal com- 
ponent of the 2nd “group is very small in 
amplitude compared with those of the others. 
This fact corresponds well to the seismic ray 
paths considered just now. 


5 


Seismic ray paths in the strata. 


iss 19) 


Emergent angles of the ray paths in Fig. 
19 can be calculated as follows; : 


Wave Groups Generated by a Very Small Explosion. 


Pty =sin-(vp;/Cory) imaginary, 


Oy1= SiN-(Vp/Cor1) 40°, 
9y1= sin! (Vp3/Cor) 10, 
9y12= Sin“(Vp»/Coq1) 25°. 


Imaginary value in the above relation seems 
to be difficult to explain at present and the 
discussion on it will be left for future. 

In the case of triple layers, unfortunately, 
few dispersion curves have been calculated as 
to M-waves. Any law as that in the case of 
doubly stratified layers can not be proved 
theoretically. Now let us assume the follow- 
ing relation, similar to the previous one, 

onal (i, A, 
4 alge Ga) 
(m means any positive integer). 


Putting the following values (see Fig. 13 and 
Table 4) 


Om =—1s3.X 10? m/s, vps=2:5 x10? m/s, 
H,=1.0m, H,=3.0m, T)=2.8x 10-? sec 
into the above equation, we get 


an—1. A e ee, 


en OA A 4 


Since the values above used might be some- 
what, though not greatly, uncertain, (227—1)/4 
will preferably be taken as 3/4 instead of 
2.9/4. Thus we have 3/4 wave length law 
also in triple layers. 

If we consider this result alone, the 2nd 
wave group seems to be ;M@,; in Fig. 17 (a). 
But we have a contradiction with respect to 
its phase velocity, because this is much larger 
than the phase velocity of »M2; which will be 
found later to spread over upper two layers. 
The theoretical phase velocity corresponding 
to the first minimum group velocity must be 
always smaller than that corresponding to the 
second one. Thus our 2nd group cannot be 
.M,,, although it satisfies the 3/4 wave length 
law. 

We have also examined the Airy phases in 
the dispersion curves calculated by Totstoy 
and Eucene (1953), and have got the relations 
given in Table 6. Although these are con- 
cerned with doubly stratified layers, they 
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Table 6. Various wave lengths law predicted 
by theoretical dispersion curves. 


To(Hjog) |) Shranch and phase 
1/4 | Mi, | 2M, | 
3/4 | My | 1M 2M)» 
5/4 | | My», 


| 


may be applied to triple layers by the same 
reasoning as above. We see M,, and ,My. 
also satisfy the 3/4 wave length law like 
M,). 

Now phase velocities of both Mj, and .Mop, 
the latter in particular, become much larger 
than that of ,M@.;. Thus we can attribute the 
2nd wave group to the phase corresponding to 
the second minimum group velocity of Mb- 
wave in the superficial two layers. 

In short, all predominant phases in our 
seismic records Fig. 1 and Fig. 2 have been 
explained, although the blast used was re- 
stricted to a cap. 


$9. Wave Groups Generated by a Larger 
Explosion 


In fact, the seismic records by 30 gm dyna- 
mite have more complex features than in the 
previous cases, the other conditions being kept 
equal to each other. Some of the records in 
the later cases are illustrated in Fig. 20, and 
the time-distance curves of every trough are 
shown in Fig. 21 where chain and broken 
lines are respectively group and phase veloci- 
ties previously observed. In both figures, 
predominant phases are marked by I’, U’,---. 
Among them III’ wave group has the most 
complex feature, and this complexity may be, 
as seen from Fig. 21, due to interferences of 
the 3rd and the 4th wave groups already 
mentioned. Variations of periods or wave 
forms with charge amounts will be discussed 
in another opportunity. In this section, we 
shall be concerned only with the classifications 
of wave groups caused by an explosion. 

The time-distance curves of I’ and II’ groups 
are just the same as those of the Ist d) and 
the 2nd (II) groups. The former groups must 
be identical respectively to the latter ones. 

Wave groups which have not been observed 
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Fig. 20 (a) U, 20.25m ~ 23.25m 
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Fig. 20 (b) #, 10.75m ~ 27.25 


Fig. 20. An example of the seismic records by 30gm dynamite. 


previously are IV’ and V’ groups marked in Admittedly the law is not satisfied perfectly, 
Fig. 20 and Fig. 21. Several data on these but their phase and group velocities seem 
groups are given in Table 7. Now assuming to show that they are the phases near the 
Arry phases of ».M., and M,,-waves in the 
Table 7. Collected data on IV’ and V’ groups. superficial two layers. Thus we have the 


wave eid | phase | wave | group results in Table 8. 
group P | velocity length | 
| got Youle 


| Pelle Pelee ce velonity Table 8. Classifications of wave 
Tyee ee 0 | 6m | 0.3.0-35 — —— nae | 
ieee Lute | e. wave Be 4 | i 
v' | 5~6 | 0.7~0.9 | 5m | 0.2~0.25 eolp 4 I maa | Ill | IV M+ TV Iv’ Vv’ 
time unit 10-2sec; velocity unit 102 m/s 
. type | P| Mx» »My, Mir | 2»Mo,, | 2M | My 
again the (272—1)/4 wave length law, we have | | | | My, | 
for the both groups layer 12 jk ld joe std 2 eke 
2n—1_ ae: 3.0 \e 1.3 Fig. 22 shows particle orbits of all phases 
4 ai Nealleca): 2els; ae observed. The reductions from the recorded 
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Fig. 21. The time-distance curves 
in the case of 30gm dynamite. 
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Fig. 22. Particle orbits of all phases observed. 


orbits to the true orbits are not made in Fig. 
22. Particle orbits of IV’ and V’ groups 
seem to contradict theoretically expected 
ones in doubly stratified layers (SmzaAwa and 
Kana (1940)). But the hodograph directly got 
from the seismic record may have some dif- 
ferent feature from the true one. We shall 
not adhere to the problem at present. 


§10. Shielding Effects of Upper Layers 
Looking at the initial motion in our records 

once more, we can find very poor waves be- 

fore the wave of velocity 18x10?m/s, which 
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is taken as the initial motion in the present 
discussion. The velocity of that very feeble 
initial motion seems to be 25x 102 m/s, and 
this may be supposed to be some refracted 
wave due to a lower stratum than that we 
have treated by now. But its amplitude is 
too small to be observed on all records. 
Looking again at the boring data in Fig. 13, 
we see that the velocity ratio of the light 
clay to the following one is only 25/18=1.4, 
whereas: that of the dark clay to the light 
clay is 18/2.5=7.2. Thus the dark clay layer 
has shielding effects for seismic waves, not 
only for refracted waves but also for surface 
waves. Therefore all wave groups observed 
by us are concerned only with strata above 
the light clay layer under which there might 
be many layers actually. By this reason, we 
have only a finite number of wave groups 
at the present experiments. 

It will, however, remain questionable why 
iM, M,. and M2, in the uppermost layer and 
VM, and M,, in the superficial two layers 
have not been observed in our seismic records. 
According to our field experiences, JM, seems 
to have smaller amplitude than ,/%, whose 
group velocity is smaller than that of :My. 
M,, and M,, in the uppermost layer must 
have periods of about 1.0 10-* sec, where the 
magnification of our recorders is very small. 
The absence of M,. in the upper two layers 
cannot be explained at present. It will not be 
explained until we have theoretical knowledges 
on relative amplitudes of M-waves in these 


layers. 


§11. Conclusion 


The possibilities of existences of many 
waves, direct, refracted, reflected and various 
surface waves, have already been predicted by 
many theoretical investigations. But it is a 
different problem to ascertain these waves on 
actual seismic records. It is most important, 
in practice, to recognize the characters of 
the waves predominant on our seismograms. 
Some of the waves are useful and others may 
be obstructive for some purposes. Theory and 
practice in general have not yet stepped up to 
each other, though a few attempts were made. 
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The present paper, the author believes, has 
given the first thorough discussions in this 
direction. 

The present model seismology is confined in 
a very small scale, but it will be extended to 
larger scale observations. The refraction 
method alone cannot decide a unique subsur- 
face model. In such a case, it may be useful 
to count even the number of the surface 
waves alone. In order to obtain a good reflec- 
tion record, the shot depth must be chosen 
lower than a shielding layer thus determined. 

The present discussions on amplitude dis- 
tributions in the strata concerned and the 
quarter wave length law are incomplete, but 


they will. be made perfect in another 
opportunity. 
Acknowledgements 


The present author expresses his hearty 
thanks to the members of the Seismic Explora- 
tion Group of Japan headed by Prof. Kenzo 
Sassa for the kind cooperations in the present 
experiments. Many thanks are also due to 
the technical experts of Teikoku Oil Company, 
Transportation Technical Research Institute of 
Japan and Geological Survey of Japan. The 
author wishes again to express his thanks to 
the Government Authority of the Ministry of 
Education for the grant of the Science Re- 


Kyozi TAZIME 


search Fund by the aid of which the present 
investigation was made possible. 


References 


HATAKEYAMA, T. 

1954 Hansha Jikken. (in Japanese) Rep. Seis. 

Exp. Group of Japan. No. 1, 3. 
KANAI, K. 

1951 On the Group Velocity of Dispersive Sur- 
face Waves. Bull. Earthq. Res. Inst., 29, 
49. 

MURAUCHI, S. 

1955 S ha o hassei suru Sochi ni tsuite. (in 
Japanese) Rep. Seis. Exp. Group of Japan, 
No. 4, 21. 

SEZAWA, K. 

1927 Dispersion of Elastic Waves Propagated 
on the surface of Stratified Bodies and an 
Curved Surface. Bull. Earthq. Res. Inst., 
3, or 

SEZAWA, K. and KANAI, K. 

1940 Dispersive RAYLEIGH Waves of Positive or 
Negative Orbital Motion, and Allied Pro- 
blems. Bull. Earthq. Res. Inst., 18, 1. 

SUZUKI, Z. and Sima, H. 

1954 On Forms of Seismic Waves Generated by 
Explosion 1. Sci. Rep. Tohoku Univ. 
Geophys., 6, 85. 

TOLSTOY, I. and EUGENE, U. 

1953 Dispersive Properties of Stratified Elastic 

and Liquid Media. Geophysics, 18, 844. 


